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ABSTRACT 

Three computer programs were written to find the eigenvalues 
and eigenfunctions of acoustic normal modes in the ocean. The 
programs used two different methods: an iterative finite difference 

scheme, and a method based upon the WKB approximation of quantum 
mechanics. The methods assume a flat fluid bottom and are designed 
for any arbitrary sound speed profile. While the results of both the 
finite difference and the WKB methods agreed, the WKB method 
proved faster. 
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L INTRODUCTION 



Since the Second World War considerable interest has developed, 
in both naval and scientific communities, in the prediction of sound 
propagation within the ocean. For the scientist, propagation informa- 
tion can be a vital part in investigating the acoustic, physical and 
sometimes chemical properties of the ocean. For the naval commu- 
nity, propagation is a vital element in the prediction of acoustic sen- 
sor performance. The prediction of such performance is not limited 
to the design and development of sensor systems, but is increasingly 
important in the operational employment of such sensors and selec- 
tion of the most fruitful tactics. 

Accordingly, there has been considerable development within the 
naval community of numerical acoustic propagation models. These 
models have, for the most part, been based upon the theory of ray 
acoustics. The techniques employed have developed considerable 
sophistication in order to deal with ray limitations with caustics, 
frequency effects, and interference. However, for prediction of prop- 
agation over great ranges and at low frequencies, ray techniques 
have two serious limitations: 

(1) The long ranges involved require long computation times. 

Each ray must be ’’traced” over many successive short time steps to 
simulate travel of the wavefront. Usually more than one hundred such 
rays are traced. Thus, as ranges increase, the computation time 
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increases proportionately. Such long computation times make the 
routine operational use of ray trace programs covering paths of 
thousands of kilometers too expensive to be practical. 

(2) Ray trace programs largely ignore such frequency dependent 
effects as diffraction and interference. The assumptions of ray a- 
coustics require that the wavelength of the acoustic signal be short 
in comparison to the depth and velocity gradient scale. As lower 
frequency signals are considered, the wavelength becomes an appre- 
ciable fraction of the depth scale, and so violates a basic validity 
assumption for ray acoustics. 

Because of these limitations there has been increased interest 
in normal mode theory as a basis for long range propagation models. 
The theory of normal mode propagation is not new and has had a num- 
ber of shallow water applications [Officer 1958, Bucker and Morris 
1965 ]. A. O. Williams (1970) has outlined a method whereby the 
normal mode technique could be applied to deep ocean acoustics. As 
a part of any such technique, both the characteristic horizontal wave- 
number (eigenvalue) and the mode shape (eigenfunction) must be 
solved for each mode. To do this, three basic methods have been 
used: 

(1) The wave equation has been solved in closed form [Tolstoy 
and Clay 1966, Williams and Horne 1967]. In order to do this, the 
sound velocity profile must be fitted to an assumed analytic function. 
Such functions have taken the form of Epstein profiles [Bucker and 
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Morris 1967], profiles with a constant gradient of the reciprocal of 
the sound velocity or sound velocity squared [Williams and Horne 
1967 , Tolstoy and Clay 1966]. The requirement that the sound veloc- 
ity profile be fitted to an arbitrary function places severe limitations 
on the flexibility of any model. 

(2) The wave equation has been vertically integrated while meet- 
ing appropriate boundary conditions. This method has the advantage 
of being adaptable to any arbitrary sound velocity profile. Kanabis 
(1972) and Newman and Ingenito (1972) have developed two such shal- 
low water normal mode models. The first program presented with 
this thesis, NORMOl, is based upon these two shallow water models. 

(3) The Wentzel-Kramers -Brillouin (WKB) approximation of 
quantum mechanics provides an analytic solution to the wave equation. 
This solution is singular at depths which correspond to ray vertex 
depths. However, the WKB approximation may lend itself to solu- 
tion of the eigenvalue, after which the eigenfunction may be solved 

by an integration similar to the previous method. The last two pro- 
grams presented with this thesis, N0RM04 and N0RM05, use this 
method. 

The long term goal of this research is to develop a deep ocean 
normal mode acoustic propagation model of sufficient computational 
speed to be considered for operational use. In order to achieve this 
goal, a rapid method of solving for the eigenfunction must first be 
developed. Based upon this need the immediate objectives of this 
thesis are: 
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(1) To compare the results of the programs employing methods 

(2) and (3) above with analytic solutions, the published results of 
other programs, and each other. 

(2) Determine whether the WKB method offers promise of im- 
provement in terms of time versus accuracy. 

(3) To outline future improvements to either method based upon 
the preliminary results. 
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II. THEORY 



A. WAVE THEORY 

In this section we apply a mathematical separation of variables 
to the wave equation; then we impose boundary conditions on the re- 
sulting separated vertical equation. The results of this mathematical 
procedure offer us a differential equation and boundary conditions 
capable of numerical solution. We will close this section with a dis- 
cussion of some general properties of such a solution. 

The following discussion is based upon the treatments given by 
A. O. Williams (1970), and I. Tolstoy and C. S. Clay (1966). For a 
more complete descriptive treatment the reader is referred to the 
discussion given by Williams. 

1 . The Wave Equation 



The simple scalar wave equation for small amplitude waves 



is 




( 1 ) 



where $ is the displacement potential. Here the displacement 



of a fluid particle from its rest position, is represented by 




( 2 ) 
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and the acoustic pressure, p, is given by 



P = -D^ = . 



( 3 ) 



Thus, if we can express ^ as a function of space, we can then find 
the resulting acoustic pressure and sound propagation loss due to 
spreading and refraction. 

If we use a cylindrical coordinate system and assume ^ is 
a function of range, depth, and time, we can rewrite (1) as 



ir / 



, 1 



(4) 



Let us further assume that $ is separable in terms of r, z, and t, 
such that 

^(r,z,t) = R(t-) Z(Z) (5) 



where co is the source angular frequency, Z{z) is the vertical dis- 
placement potential function, and R(r) is the radial displacement 
potential function. We immediately notice that 



= - CO" $ . 

Using equations (5) and (6) we rewrite equation (4) as 



( 6 ) 



^ ^ f -p . 4 b" Z _ Q 

r R arl ar / H d 



(7) 
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This expression has a range dependent term, a depth depen- 



dent term and a term with c, the sound velocity. If we could assume 
that c is independent of either depth or range, the expression is sep- 
arable into depth and range differential equations. The obvious choice 
is to assume that c varies only with depth. This assumption, termed 
horizontal stratification, although not true over great distance, is 
commonly made in oceanography and ocean acoustics. Actually, for 
our purposes it will be sufficient if c is horizontally stratified in the 
vicinity of our solution, and further if any horizontal variation of c is 
small with respect to the vertical variation. Making such an assump- 
tion, we continue by separating (7) into range and depth differential 
equations 



Solutions to equation (8) include the Hankel functions of first 
and second kind representing cylindrical waves. 



For distances such that kj.r is much greater than one, this function 
can be approximated by its asymptotic form 




( 8 ) 




( 9 ) 




( 10 ) 




( 11 ) 
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We now have an assumed time dependence and a solution to the range 
dependent function of ^ . Thus, at long range from the source we 
have for $ 






(z) 






•I 



4 



( 12 ) 



2 

Here kj. , an eigenvalue which appeared in the separation process as 

a mathematical constant, will be seen to be a horizontal wave number. 

We will consider only positive values of k^., which represent outgoing 

wavefronts. It remains for us to find a method of solving for Z(z) in 

2 

terms of the separation parameter, k^ . 

We can rewrite (9) as 









(13) 



This equation is the separated, space form of the wave equation, or 
Helmholtz equation, for the vertical wave component. The expression 
in parentheses represents a vertical wave number, k^, and is related 
to the (true) wave number, k, and horizontal wave number, k^, by 



cZ 



(14) 



Equation (13) with a pair of associated homogeneous boundary condi- 
tions forms a Sturm -Liouville problem, which is solved in terms of a 
set of eigenvalues and corresponding eigenfunctions . It now remains 
for us to define and employ boundary conditions associated with the 
ocean vertical sound profile. Before so doing, however, we will di- 
gress a moment to introduce a notational convenience. 
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■ft. 



Since the value of c does not vary greatly (perhaps five per- 
cent) throughout the depth profile, the square of the vertical wave 
number represents a small difference between two nearly equal num- 
bers 






(15) 



In order to conserve accuracy and facilitate computational speed, a 

notational convention is borrowed from quantum mechanics. The 

2 

square of the (true) wave number, k , is subtracted from an arbitrary 
constant of the same order of magnitude (here the maximum possible 
wave number) to form a ’’potential function. ” 



= # • 

' ' ^ 2 . -max r^. 



C ^ 



— 

-max 



(16) 



The square of the horizontal wave number (our separation constant) 
is also subtracted from the maximum wave number to form a quantity 
corresponding to the ’’energy level” of quantum mechanics, and our 
new eigenvalue 






( 17 ) 



With this notation, the Helmholtz equation (13) becomes 

-[e-y<.)]z=o , <'8> 

a form similar to Schroedinger *s equation. By this notation we hope 
to facilitate the adaptation of the results of quantum mechanics to the 
problem at hand. 
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2. Boundary Conditions 



We continue by developing appropriate vertical boundary 
conditions . 

Figure 1 represents a typical deep ocean sound profile. 
Figure 2 is the same profile represented in terms of V(z). Note that 
at some intermediate depth the value of sound velocity reaches a min- 
imum at the sound channel axis. Note also that the function V(z) 
forms a ’'potential well” with depth, representing the deep sound 
channel. Our domain of interest is bounded by the sea surface and a 
flat bottom. 

The air -water boundary at the surface approximates a free 
surface at which the stresses vanish. In a fluid this boundary corre" 
sponds to 



Z(o) = 0. (19) 

If the limit of the second derivative of the vertical displacement func- 
tion is finite at the surface, the stipulation that V(z) dis continuously 
approaches infinity at the surface has the effect of making Z(z) vanish 
at the surface. Thus, in terms of our quantum mechanics analogy, 
the free surface boundary condition corresponds to a perfectly rigid 
wall of the potential well. 

The ocean bottom provides a boundary more difficult to 
specify. The ocean bottom has a complex, often layered, horizontal- 
ly varying structure. In addition, little is known of the specific 
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TYPICAL DEEP OCEAN SOUND PROFILE 
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TYPICAL DEEP OCEAN VELOCITY FUNCTION 




Figure 2 
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acoustic properties of much of the bottom required tomodel fully the 
bottom boundary conditions. In the face of such complexities we will 
approximate the bottom with a relatively simple model. 

We will assume that the bottom is represented by a homo- 
geneous fluid of infinite depth, and constant velocity. Across a fluid- 
fluid interface we require that both pressure and vertical particle 
motion be continuous. From (3) and (6), continuity of pressure re- 
quires 

e,Zo = (20) 



at the bottom. From (2) we find that continuity of vertical particle 
motion or displacement at the bottom requires 



^Zo ^ 
d z 



( 21 ) 



We combine equations (20) and(21) to form a single bottom boundary 
condition 



1 J_ 

fo Zo PtZw 



( 22 ) 



Since we have specified the bottom to be of constant velocity 
and density, a known solution to equation (18) for the bottom region 
is 

Hb= ^ (23) 

where Kg, a real constant, is the imaginary bottom wave number, 
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Ks-J V, - E' . 



(24) 



If the amount of energy represented by the bottom solution is to be 
finite, the limit of the integral of from the bottom boundary to 
infinity must in turn be finite. This condition requires that the expo- 
nential growth term of be suppressed. Thus, (23) becomes 

•"7 -KqZ 

Zb = Cae , (25) 



and (22) becomes 



Jl = 
Zo ^2 




( 26 ) 



We have now specified two boundary conditions which will 
enable us to find particular solutions to the vertical differential wave 
equation (18). Before proceeding with the technique of finding such 
solutions, let us discuss some of their properties. 

3. The Nature of the Solution 

Before we find the solutions, Z(z), to the differential equa- 
tion (18) and boundary conditions (equations 19 and 26), it will be 
beneficial to consider the form and properties of the solutions we 
expect. 

a. In the previous section we saw that the solution in the 
bottom region is of the form of a decaying exponential (equation 25). 
This solution requires that the quantity Kg in equation 2 5 must be 
real. This requirement in turn implies that 
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E<Vb 



(27) 



or 



c 



(28) 



In order to have propagating waves, it must also hold that 



E>0, 



(29) 



or 



kj.< 3^ 



(30) 



min 



Thus, we have an upper and lower bound for E and k^. 



V^>E>0 



(31) 



and 



C , C 

min b 



(32) 



b. Within these limits the boundary value problem is so con- 
strained that it can only be solved in terms of a finite number of 
eigenvalues, and corresponding functions E^(z). This eigen- 

function, properly normalized, forms what is termed a normal mode. 
At sufficient range a vertical sound pressure profile can be approxi- 
mated by a linear combination of such eigenfunctions. 



NJ 

pCz) = Z] Cm . 



(33) 
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TYPICAL MODE SOLUTION 




Figure 3 
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c. At depths such that V(z) the eigenfunction Z^(z) is 

oscillatory in character. We can represent such a function by 



Z(’2) ^ Su\(S(£>) 



(34) 



where ^(z) and S (z) are undefined functions of depth. We shall 
refer to such depths as the oscillatory or ’’real” region (a reference 
to the fact that the vertical wave number, k^, is a real number). At 
depths such that E^< V(z) the eigenfunction Zj^(z) becomes quasi- 
exponential in character, with a form we can represent by 



where, again, ^ (z) and S(z) are undefined functions. We shall refer 
to such depths as the exponential or '‘imaginary region. 

The depth at which the sign of [e^-V ( z^ changes is termed 
a turning point. At the turning point the character of the eigen- 
function changes from oscillatory to quasi -exponential (Fig. 3). 

d. Each mode, corresponds to one or more rays which traces 
paths within the oscillatory region between turning points, surface, 
or bottom boundaries. The angle the ray makes with the horizontal 
is defined by 



Cos = CCZ) 




(36) 



From this we see that the mode turning point corresponds to the ray 
vertex depth (that depth at which a ray reaches its maximum depth 
excursion). When the oscillatory region of a mode is bounded by the 
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free surface boundary or the bottom, the corresponding ray is sur- 
face or bottom reflected, respectively. 

e. The lower boundary placed on kj. represented by equation (28) 
has a more familiar and perhaps more satisfying physical meaning. 
Substituting equation (36) into (28) we have 



Cos (A > 

^ Cb 

At the bottom this is the expression for the critical angle. Thus, we 
are limiting outselves to consider only those modes whose corres- 
ponding ray strikes the bottom at a grazing angle less than the criti- 
cal angle. Such modes are widely called ''unattenuated modes." 

Also, there exists an infinite set of modes such that 






(38) 



However, because of bottom reflection loss, those modes tend to 
attenuate rapidly with range. For propagation problems at a con- 
siderable range, the "attenuated modes" contribute an insignificant 
amount to the acoustic pressure field, and are ignored. 

f. Given the oscillatory nature of the eigenfunction Z^, we see 
that the eigenfunctions for each successive mode must change sign 
one additional time (Figure 4). Each sign change will be referred 
to as a mode crossing. Thus, corresponding to each eigenvalue E^, 
there corresponds an eigenfunction with m-1 mode crossings. 
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Figure 4 
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4. Mode Parameters 



We should now introduce two sets of parameters, derived 
elsewhere, which are used in normal mode calculations, 
a. Phase and Group Speed 

Phase speed describes the horizontal speed of the wave 
front represented by a mode. As the speed of advance of a wave- 
front of constant phase, the phase speed is given by 



= f ■ ( 39 ) 

Group speed is the rate of energy transport in the hori- 
zontal direction. Tolstoy and Clay (1966) give an expression for 
group speed based upon a theorem by Biot (1957). This method cal- 
culates group speed by 



Cg- 1 -i- 

Lp O' 



(40) 



where V' is the normalization integral, 




and 



(41) 




(42) 



b. Bottom Decay Coefficient 

Kornhauser and Raney (1955) give an expression which 

calculates the effect of a bottom absorption on the mode amplitude. 

Assume the wave number in the bottom is complex and given by 

30 



( 43 ) 



where ^ is a bottom attenuation coefficient. Then the horizontal 
wave number must also be complex and given by 

Hr - i-c'v • (44) 

If B is small with respect to ^ , then the ratio of c< to (!? for a given 

^b 

mode is approximated by 

( 45 ) 

bolbe>m 

B. THE WENTZEL-KRAMERS-BRILLOUIN (WKB) APPROXIMATION 
The WKB approximation is a method borrowed from quantum 
mechanics which will enable us to approximate the eigenvalues, ^^9 
of equation (18). In the WKB approximation we assume a form for 
the mode solution which is valid at depths removed from the turning 
points. We then develop a characteristic equation for the eigenvalue, 
based upon the assumed solution. This section follows in many re- 
spects the descriptions given by Tolstoy and Clay (1966) and Lauvstad 
(1971). 

1 . The Characteristic Equation 

Based upon the sign of [E - V(z)] the vertical wave equation 

(18) applies to two domains, which we have termed the 'imaginary” 

region and *'real” region. Let us rephrase equation (18) into two 

separate equations for the two domains. So doing, we have 
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for E y V(z), 



(46) 



zi' * Z. ■ 0 

and 

Zz - K/ Zz =0 for E < Vfz). (47) 

Here we have defined 

" [E * V( 2 )] when E y V(2), (48) 

and 

Kj = [Vc 2 ) • E] when E < Vc2). (49) 

Now let us assume a form for the solution to equation (46), 
similar to that which we have previously postulated in equation (34). 
The assumed solution is 

. = <50) 

By substituting equation (50) into (46) we have as a real part 

S'«) - 5«) ) = 0, ^ ) 

and an imaginary part 

v5(*)S''?z) +25‘^(z)S'z) = 0. (52) 

This equation immediately yields an expression for § (z), 

5fz) - A/S^ . (53) 
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Now let us simplify equation (51) with an assumption. Assume that 
the first term of equation (51) is negligible with respect to the others, 
or 



1 



6 (£) 



«1 . 



(54) 



With this assumption we obtain as an expression for S(z) 



S'f.," . Jk} , 



(55) 



or upon integrating 






(E) = - j Jit dz 



(56) 



We now have as the assumed solution 



Zi(H)= |VCos Si(z) + bSin 



(57) 



or, in a more convenient form, 



HifH) * a Sin + ^ 



(58) 



Similarly we obtain as a corresponding solution for the 
"imaginary" region, 



2,<z. 



(59) 



where , 82 ( 2 ) is defined by 



SgfH) ~y Kj dz. , 



( 60 ) 
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Having written the two solutions, let us take a closer look 
at the assumption implied by equation (54). With equations (53) and 
(55) we can rewrite (54) for the "real" region as 



1 




d lod 

c/2 






«l , 



( 61 ) 



and for the "imaginary" region, 






dz > 



( 62 ) 



This assumption requires that the order of magnitude of the vertical 
wave number, and thus the sound speed, vary slowly with respect to 
depth. This is generally the case in the ocean, where the total sound 
speed variation with depth is on the order of five percent. We also 
see that our solutions lose any validity as z approaches a turning 
point. Recalling the analogy of the turning point to the ray vertex 
depth, we see that as a wave approaches the turning point its orienta- 
tion becomes horizontal. Near a turning point the vertical wavenum- 
bers, k^ and approach zero and the conditions of equations (61) 

and (62) are no longer satisfied. Thus the solutions represented by 
equations (58) and (59) become singular at the turning point location. 

The WKB method does not, at this point, appear satisfactory 
for the computation of the mode profile or eigenfiinction How- 

ever, if the WKB approximation can be used to obtain an accurate 
estimate of the eigenvalues, a finite difference scheme can then 

be used to calculate the eigenfunction, Z^^. 
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Consider a typical V(z) profile, with upper and lower turning 
points, as well as bottom and pressure release boundaries (Figure 5). 
Remember that we have specified a derivative boundary condition at 
the bottom. This boundary condition determines the eigenfunction 
(except for a constant of multiplication) in the "imaginary'' region 
ZTLT^Zb* lower turning point, the ''imaginary" region 

solution, Z 2 (z), corresponds to a particular "real" region solution, 
Zj(z). We can designate such a solution by the use of an initial phase 
angle, 9 such that 

_i_(C + De — *-A. Sin fSi(£)+ 0u) (63) 

/k; » v ; 

where 

= /KjJh , (M) 

tfoTTDM 

and 

r z 

Sjcz) = J . (65) 

In a similar manner we can stipulate that in order to satisfy the free 
surface boundary condition, Z(0) = 0, the "real" region solution, Zj^ 
must correspond to a particular "imaginary" region solution, Z 2 . 
Again we can represent this as a phase angle, 0g, in the "real" region 
solution, Zj. Thus we have required the argument of Sin (S ^ (z) 
at the upper turning point (or surface) to be our specified value 
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( 66 ) 




If we let 0u ** 9s rewrite this equation as 




Since this equation describes the arguments for which the sine of the 



integrated vertical wave number is zero (a Dirichlet boundary condi- 
tion), the surface boundary condition will also be satisfied if 



^TL 

Values of E satisfying this characteristic equation will be the eigen- 
values, for which we are searching. Note that this equation is 

also that given by Tolstoy and Clay (1966) as the characteristic equa- 
tion for stratified acoustic waveguides. 

Now that we have a characteristic equation, some method of 
evaluating and 0^ is required. The WKB approximation provides 
a method for finding the eigenvalues, E^, once 9^ and 9^ are evalu- 
ated. 

2. Turning Point Connection 

The most vexing problem in using the WKB approximation is 
connecting the two solutions (58) and (59) across the turning points. A 
number of techniques have been developed, and here we will consider 
two such techniques which we will term the asymptotic and Airy phase 




111 = 1 , 2 , ( 68 ) 



methods . 
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a. Asymptotic Method 

We have two corresponding solutions represented by 
equations (58) and (5 9) which correspond to the same function on either 
side of a turning point. As the turning point is approached let us re- 
quire that the limit of the first derivative divided by the mode value 
be continuous. This is in fact the same as the fluid -fluid boundary 
condition of equation (22). For the purpose of these turning point 
discussions we will designate the turning point depth, equal to 

zero and corresponding to Zq of equations (56) and (60). Further, 
z will be positive towards the ’’real" region and negative towards the 
"imaginary’’ region (Figure 6). Thus, from equation (22) we have 



ixm E' = Z' 

e-o^.2i 2-0- 



(69) 



By taking advantage of the approximation involved in equation (61), we 
get by substitution for the above 



Lm Jk, Cot(Si(E)-0o) = hyn Kh 



Taking the limit we have 

k = CotCOo) = D - C ■ (71) 

D * C 

It should be noted that this approach is simpler than the Airy phase 
technique (described next) and may not be as accurate. It is included, 
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however, since it provides an interesting comparison with the Airy 
phase method in the computer programs, 
b. The Airy Phase Method 

What follows is an outline of a technique developed by 
V. R. Lauvstad (1971) of the SACLANT ASW Research Centre. 

Lauvstad assumes asymptotic solutions of the same 
form as equations (57) and (59) in regions remote from the turning 
points. The gradient of the square of the vertical wave number a- 
cross the turning point is assumed to be linear, such that 

( 72 ) 

This reduces the vertical Helmholtz equation to 

One solution to this equation is the Airy phase function 



Z = AAt (-S(z)) BBx(-Scz)), 



(74) 



where S(z), a depth integration function, is positive in the ’’real" 
region, and negative in the "imaginary" region. The asymptotic 
forms of the Airy phase are then compared with our corresponding 
assumed solutions. We can represent the asymptotic form of the 
Airy phase for negative real argument (the "real" region) as 

Sin S(z) + (A^B) Cos 5(2)] , (75) 
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and for positive real argument as 




1 [a 
p^iz 




( 76 ) 



By matching the above coefficients with those in our solutions we have 
a set of coefficient matching equations representing the required turn- 
ing point connection. 



a = + C) 

n? 



h-- i(ED-C), 

/T’ 



( 77 ) 



C=1 (CL-b) 

W 



D = l(a4- b) 

nr 



( 78 ) 
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III. THE PROGRAM METHODS 



Here we will discuss the basic methods used by the three pro- 
grams (NORMOl, N0RM04, and N0RM05) presented as a part of 
this thesis. This section is intended as only a brief introduction to 
the programs and their characteristics. A following section will give 
an outline of the program algorithms. 

A. NORMOl - A FINITE DIFFERENCE TECHNIQUE 

NORMOl is essentially an adaptation of two shallow water pro- 
grams developed by Kanabis (1972) and Newman and Ingenito (1972). 
For a given value of the horizontal wave number, k-,.^or E, the pro- 
gram begins with the bottom boundary condition (equation 22), and 
steps through a finite difference scheme to the surface. The finite 
difference scheme, which is derived in appendix A, is from Kanabis 
(1972). The finite difference scheme is a third-order Newton forward- 
difference scheme which iterates both the eigenfunction and its deriv- 
ative. Using this difference scheme, a search is made for those 
values of k^, or E, which most closely satisfy the surface boundary 
condition (equation 19). Such values are the sought eigenvalues. 

A shortcoming of this technique is that under certain circum- 
stances the eigenfunction at the surface, tends to grow expo- 

nentially rather than to decay towards zero. Some mention of this 
phenomenon is made by both Kanabis (1972) and Newman and Igenito 
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(1972). Kanabis recommends decreasing the incremental search 
limits for the eigenvalue (vhich imposes a consequential increase in 
computation time) to overcome the phenomenon. Newman and Igenito 
set the mode value equal to zero at depths above the last mode cross- 
ing or minimum. It happens that the circumstances for this expo- 
nential growth, or degenerate solution, occur frequently in deep 
ocean acoustic profiles. Thus, a more detailed treatment of the 
phenomenon is required. 

To understand this degenerate solution let us consider a typical 
deep ocean sound profile and our general forms for the eigenfunction 
(equations 34 and 35). Consider a mode whose eigenvalue, E^, inter- 
sects the velocity function, V(z), at some depth, below the sur- 

face. This is a common (in fact, the usual) occurrence for the lower 
modes in the deep ocean, where a deep sound channel and thermocline 
usually exist. As noted by Schiff (1955) and Lauvstad (1971), if the 

value of S(z) in equation (34) is not exactly^ at the turning point, the 

4 

coefficients C and D of equation (35) are both non-zero. Thus, the 
exponential growth term has some finite value. Given both sufficient 
depth between the turning point and the surface, plus sufficiently 
large values of [V(z)-E], the growth exponential will eventually be - 
come dominant. Thus, the solution tends to grow exponentially as the 
surface is approached. 

This degenerate solution can complicate the search for eigen- 
values if it is based solely upon the value of Z(0). With the 
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degenerate solution, the value of Z(0) will oscillate between extreme- 
ly large positive and negative numbers, even with only small changes 
in the horizontal wave number or E. For the lowest modes of a deep 

sound channel profile the values of Z(0) can so oscillate with incre- 

- 14 

mental changes in E of only one part in 10 . Thus, a method such 

regula falsi (described later) can become inappropriate since a 
map of Z(0) as a function of E can appear discontinuous (Figure 7). 

In order to circumvent this difficulty, NORMOl uses a halving 
procedure based upon the number of mode crossings. When the sur- 
face boundary condition is met, the eigenvalue E^ is the largest 
possible value of E such that the eigenfunction has m-1 mode cross- 
ings. Based on this property, the number of mode crossings is used 
to converge upon the mode eigenvalue. The method of successive 
halving, although not sophisticated, has the versatility to deal suc- 
cesfully with the degenerate solution and any arbitrary profile. Once 
an eigenvalue is found, the function Z^ is then checked for a degen- 
erate solution near the surface. If it is present a new profile is 
iterated in the upper ''imaginary'* region and matched to the lower 
profile. This procedure may cause a discontinuity in the derivative 
of at the turning point; however, its ill effects would usually be 
less serious than those of the degenerate solution. 

B. N0RM04 AND N0RM05 - THE WKB APPROACH 

In both N0RM04 and N0RM05 the vertical wave number is inte- 
grated between the upper turning point (or surface) and lower turning 
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point (or bottom). The programs differ only in the turning point 
connection formulae. Both programs integrate the vertical wave 
number using the trapezoidal rule over the depth grid points. The 
eigenvalue is converged upon using the regula falsi method. 

This method makes successive estimates of the zero value of a func- 
tion by the linear interpolation of a negative and positive value. The 
method is represented by 

Eui = E- + (£^-E-){-AS-) , (79) 

(AS^ - AS-j 

where the subscripts -f and - correspond to the values associated 
with the last positive and negative values of A S; and AS, the differ- 
ence function, is defined by 

AS - J J^^dz + 0u " •m^ • (80) 

If the mode is not found after a set number of iterations, the regula 
falsi method is dropped for a halving process. 

The phase effect of the upper and lower turning points, 0^ and 
are evaluated using the asymptotic method for N0RM04 and the 
Airy phase method for N0RM05. The vertical wave number is integ- 
rated from the uppermost turning point (or surface, if surface re- 
flected) to the lowest turning point (or bottom, if bottom reflected). 
Both programs must then provide, in addition to the values of 0^ and 
0j^, the effect of any intermediate ‘’imaginary" region or barrier. 

This occurs when multiple sound channels are encountered (Figure 8). 
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Details of the derivation of the specific formulae for 0^, 9j^, and the 



phase effect of an intermediate barrier are given in appendices B and 
C for N0RM04 and N0RM05 respectively. 

If the mode is surface reflected or bottom reflected, both pro- 
grams use the same formulae for the phase effect of the reflection, 

9^ and Since the free surface requires Z(0) = 0, for the surface 

reflected wave, we have from equations (63) and (66) 



For the bottom reflected mode, we require the conditions of a fluid - 
fluid boundary to be met (equation 22). Substituting equations (63) 
and (25) into (22) we have 



where is evaluated at the water side of the bottom boundary. 

The occurrence of multiple sound channels or ducts makes the 
WKB method more complicated. If the intermediate ’’imaginary" 
region or barrier is of sufficient extent and magnitude, sound 




TL 



(81) 



and thus 



0U - 0 . 



(82) 



ji^.CotceuV- , 

fb 



(83) 



or 






(84) 
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propagation in each channel becomes independent of the other. When 
this is the case, the waves in one channel undergo an effectively total 
reflection by the barrier. Thus, complete sound ducting occurs 
entirely within that channel. Lauvstad (1971) gives the ratio of the 
amplitudes of the reflected and incident waves in a duct as 



where L is the integral of the imaginary vertical wave number across 
the barrier, 



This result can be derived using either the asymptotic or Airy phase 
connection formulae. As L increases, rapidly approaches unity, 
and near total reflection is achieved by the barrier. The WKB pro- 
grams N0RM04 and N0RM05 evaluate L, and if it is of sufficient 
magnitude, the programs treat the mode propagation in each channel 
independently. The programs accommodate a maximum of one main 
sound channel (which contains the velocity minimum), and one upper 
and one lower sound channel. 




( 85 ^ 




( 86 ) 
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IV. THE BASIC ALGORITHMS 



A. NORMOl 

1 . General Description 

NORMOl iterates the wave equation vertically from the 
bottom to the surface using a finite difference scheme over a depth 
grid of up to 501 points. The program is written in IBM FORTRAN 
IV with subroutines used as major functional blocks. This modular 
programming is designed for both ease of modification and debugging. 
A description of the sequence of events follows. 

NORMOl is basically composed of three nested loops: a 

frequency loop, a mode loop, and an iteration loop. The input data 
are read using an input subroutine. The required input data are 
listed in appendix D. Once the data are read and appropriate con- 
stants computed, the sound velocity profile is linearly interpolated 
for all grid points. Here the frequency loop begins, the appropriate 
frequency is selected, and the velocity function V(z) is computed for 
the grid points. Next the maximum and minimum horizontal wave 
numbers and maximum number of modes are computed. If the fre- 
quency is below cut-off (the maximum number of modes is zero), a 
new frequency is selected. At this point the mode loop begins, and 
the following steps are repeated for the first to the highest requested 
or possible mode. 
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First, a method of successive halving is used to find two 



values of the horizontal wave number k^. which have m~l and m mode 
crossings. The desired eigenvalue must lie between these two values. 
Then, again using a halving process, the program converges upon the 
eigenvalue k^. The mode is considered found when one of two con- 
ditions (regulated by an input parameter, lEX) is satisfied; 



|H»,| i IHU.10 



-lEX 



(87) 



or the horizontal wave number changes between iterations by an 
amount such that 



V. - K 



l-n 



K.. 



< 10 



( 88 ) 



Once the mode is found, the eigenfunction is checked for a 
degenerate solution (exponential growth near the surface). If a de- 
generate solution is present, a new function is iterated from the sur- 
face to the upper turning point, then joined with the lower function. 
The mode is then normalized to its absolute maximum value, and the 
profile within the homogenous bottom is computed. The computed 
mode is plotted using a printer plot routine. A subroutine then com- 
putes the normalization integral represented by equation (41), and 
the group speed, phase speed and mode attenuation coefficient. At 
this point the mode and frequency loops end. 

Included in NORMOl is an output subroutine which provides 
a punched card format listing of the eigenvalues, mode parameters. 
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sound velocity profile and mode profile. The output subroutine was 
written to provide local input to CALCOMP plots and propagation 
loss programs . 

2. The Subroutines 

The following is a brief description of the major subroutines 
which comprise the NORMOl program. 

INPUT 1 provides the input data for the program. A user 
supplied page heading is read on the first card, followed by the run 
parameters, list of frequencies desired, and a maximum of 50 depth 
and sound velocity pairs. If the parameter representing the number 
of depth and velocity pairs, NUMV, is negative, all length measure- 
ments read are converted from feet to meters. This subroutine also 
prints the input data for reference. 

INTRPL linearly interpolates the input depth and sound 
velocity pairs to compute the grid point values of sound velocity. If 
required, this subroutine also extrapolates the last sound velocity 
to the bottom assuming an isothermal, isohaline gradient (0. 0 17sec”^). 

MAXMIN computes and checks the maximum and minimum 
values of the horizontal wave number, and finds the number of modes 
represented by the minimum horizontal wave number (the maximum 
number of modes). If the maximum number of modes is zero, the 
frequency is below cut-off and the program selects a new frequency 
or terminates the calculations. If the maximum possible number of 
modes is below that requested, this parameter is changed and noted. 
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HALF searches for two values of the horizontal wave number 



with m (the mode being searched) and m-1 mode crossings. These 
values represent lower and upper bounds on which are then suc- 
cessively narrowed. At the same time this subroutine ’’maps” the 
values of versus mode number in order to facilitate the search 
for later modes. 

ITERAT, the heart of the program, iterates the eigenfunc- 
tion from the bottom to the surface using equations (97) and (101). 

In addition, this subroutine counts the number of mode crossings and 
finds the maximum absolute value of the eigenfunction. 

DNORMl normalizes the eigenfunction with respect to its 
maximum absolute value. 

BOTTOM creates an exponentially decaying bottom mode 
profile from the bottom to an input depth, PLTMAX. 

MODPLT and CZPLOT are printer plot subroutines which 
plot the mode shape and sound velocity profile. Both these subrou- 
tines use a Naval Postgraduate School utility printer plot subroutine, 
UTPLOT. 

FIXUP tests for a degenerate solution near the surface. If 
the degenerate solution exists, this subroutine applies the finite dif- 
ference iteration scheme starting from the surface downward to the 
upper turning point. This upper mode shape is then scaled to join 
the original lower function at the turning point. If required, a new 
value for the absolute maximum of the eigenfunction is then found. 
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INTEGR computes the mode normalization integral (equa- 



tion 41) and the integral represented by (equation 42). The phase 
velocity is computed; then the two integrals are used to compute the 
group velocity. Finally, the mode attenuation coefficient (equation 
45) is calculated. 

B. N0RM04 AND N0RM05 
1 . General Description 

N0RM04 and N0RM05 are identical with the exception of the 
turning point connection formulae (appendices B and C explain these 
differences). Similar to NORMOl, they consist of three nested loops 
a frequency loop, a mode loop, and an iterative loop represented by 
the subroutine SEARCH. In addition, the input, interpolation, plot- 
ting, maximum /minimum and mode integration functions are per- 
formed by subroutines similar to those used in NORMOl. 

The sequence of events in the WKB method programs is 
similar to NORMOl. N0RM04 and N0RM05 read the input data and 
then interpolate the sound velocity profile over a maximum of 1001 
depth grid points. If requested, a sound velocity printer plot or 
punched card format sound profile is provided. The frequency loop 
then computes the velocity function V(z) and checks for maximum 
horizontal wave number, minimum horizontal wave number, fre- 
quency cut-off, and the maximum number of modes available. At 
this point the mode loop begins and the search for the eigenvalue is 
controlled by the SEARCH subroutine. 
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After the mode is found, the mode shape is computed using 
a finite difference scheme. The finite difference scheme is fitted be- 
tween the upper and lower turning points, then a solution which decays 
away from the turning points is calculated. After the mode shape has 
been calculated and normalized, the mode is integrated to provide the 
normalization integral, phase and group speed, and the mode attenu- 
ation coefficient. Finally, if requested, a printer plot of the mode 
shape is provided. At this point the mode and frequency loops end. 

Because of their importance to these programs, the SEARCH 
and WKB subroutines deserve more extensive description. 

2. SEARCH Subroutine 

The SEARCH subroutine controls the iterative search for the 
eigenvalues k^. The subroutine first makes an initial guess of an 
upper and lower bound for the eigenvalue. The WKB subroutine is 
then called to compute the difference function given by equation (80). 

If the difference function for the lower bound is not positive, the 
bound is successively relaxed until a positive value is found. 

The next guess for the eigenvalue is calculated using a 
regula falsi zero-finding technique (given by equation 79). The re- 
sulting incremental change in eigenvalue and the absolute value of the 
difference function is checked to test whether the desired accuracy 
has been achieved. The WKB subroutine is then called to find the 
difference function for the new trial eigenvalue. If the resulting dif- 
ference function is negative the upper bound is replaced; if it is 
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positive the lower bound is replaced. A new guess for is then 
estimated, and this iterative process continues until one of two con- 
ditions is met. The mode is considered found when the difference 
function falls below an absolute limit, 



\AS\ < 10 ( 89 ) 

or the value of the horizontal wave number approaches machine 
accuracy, 




< 10 



( 90 ) 



If after a set number of iterations, the regula falsi method has failed 
to converge upon a solution, it is replaced by a cruder and slower 
(but more versatile) halving process. 

If more than one sound channel exists, the program must 
determine if the two channels are independent of one another. If the 
L of equation (86) is of sufficient magnitude, the subroutine SEARCH 
computes the eigenvalues for each channel separately. Thus, after 
an eigenvalue has been found, two tests are made. First, it is deter- 
mined whether an upper or lower (separated) sound channel exists. 
Then, if an upper or lower channel exists, it is tested to determine 
whether it will propagate an additional mode. If such a mode can be 
propagated, the SEARCH subroutine computes its next mode within 
that channel. 
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3. WKB Subroutine 



The purpose of the WKB subroutine is to integrate the verti- 
cal wave number between the upper and lower turning points. The 
subroutine consists of a loop which steps through the depth grid from 
the surface to the bottom. In any '’imaginary” region the imaginary 
vertical wave number is integrated for use with the turning point 
connection formulae. Each time the depth loop exits an ’’imaginary” 
region for a ’’real” region, either the upper turning point phase 0^ 
or the phase effect of a barrier must be calculated. If the barrier is 
considered to be separated (if L in equation 86 is larger than a set 
value) the integration is stopped with the upper channel. 

When the depth loop reaches the bottom, the bottom phase, 
0j^, is computed for either a bottom reflected or a decaying solution 
as appropriate. Finally, the difference function is calculated and 
program control returns to the calling program. 
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V. RESULTS 



All three programs were run with a variety of sound velocity 
profiles and the calculated eigenvalues examined. The profiles select- 
ed include the published results of other models and profiles with 
known analytic solutions. In all the test runs the input variable lEX 
was set at 10. Thus, the iterations were halted when the absolute 
relative value of the eigenfunction at the surface, Z(o) > the dif- 
ference function, ^ S , (equation 80) achieved a value less than 10 
Additionally, computation was halted if the eigenvalue was incremented 
less than 10 during a halving process (in other words, when the 
solution approached machine accuracy). 

A. TEST CASE ONE - NEGATIVE GRADIENT 

This test case is one given by Kanabis (1972). The profile is 
characterized by a simple negative velocity gradient in a shallow 
water channel (Figure 9). The profile has a fast fluid bottom of the 
same density as the water. The test results are computed for a fre- 
quency of 1000 Hz. Table 1 lists the eigenvalues given by Kanabis 
and computed by NORMOl, NORM04 and NORM05 in yd \ The list 
includes a set computed by the Kanabis program and a set computed 
by an unpublished program by C. L. Bartberger and L. L. Ackler. 
Except for the higher modes of NORM05, the results from all three 

5 

programs agree to the fourth decimal place or one part in 10 . 
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TEST CASE ONE SOUND PROFILE 
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TEST CASE ONE MODES 



mode 1 MODE 2 MODE 3 
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RESULTS OF TEST CASE ONE 
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Horizontal wavenumbers in yard 



Table 2 gives the NORMOl, N0RM04 and N0RM05 results to 
the accuracy computed and relative differences between NORMOl and 
the other two programs. NORMOl required 56.23 sec of central 
processing unit (CPU) time and 412 iterations to arrive at the solu- 
tions. N0RM04 and N0RM05 required 23.50 seconds with 83 itera- 
tions and 22.74 seconds with 79 iterations respectively. The CPU 
times cited include the time required for all operations including in- 
put, sound speed profile manipulation and graph printing. 

B. TEST CASE TWO - POSITIVE GRADIENT 

The second test case is a shallow water example cited by Newman 
and Ingenito (1972). The profile is characterized by a slow, variable 
positive sound speed gradient over a fast, dense bottom (Figure 11). 
The computation frequency was 700 Hz. It should be noted that the 
accuracy parameter used by Newman and Ingenito was an absolute 
value of less than 0.001 for the surface eigenfunction Z(0). This 
roughly corresponds to a value of lEX between 3 and 4 in the NORMO 
programs. With the exception of NORM05, all the programs agreed 
to at least the third decimal place (Table 3). Considering the accu- 
racy limit set for the Newman and Ingenito run, this agreement is 
about as good as may be expected. NORMOl required 27.43 seconds 
and 196 iterations to complete computations, while NORM04 and 
NORM05 required 14. 73 seconds for 77 iterations and 14. 92 seconds 
for 79 iterations, respectively. 
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NORMO SOLUTIONS RELATIVE DIFFERENCE 



o 



X 

h 

Q 



ID 

(M 



ro 



CO 



CO 



(M 



vD 



O 

(NJ 



m 

O 

O. 

z 



CO 


r- 


00 


o 


CO 


<— H 


(M 


ID 


CM 


iD 


o 


r- 


r- 


vD 


vO 


00 


o 


00 


iD 


CO 


O 


00 


o 




(M 




CO 




00 


nD 




O 


CO 


1— H 


00 


r-- 


ID 




CO 


(— H 


o 




00 


vD 



00 

r-- 



CO 



r^ 

r^ 



CO 



vO 



CO 



nO 

o 



X 



ID 

(M 



LD 



CO 



CO 



CO 



CO 



DJ 



ID 

O 



CO 

CO 



r^ 

(M 



O 

05 

O 



CO 




00 


o 


CO 


o 


1— H 


00 


00 


vO 


vD 


r-- 


r-- 


vD 


vO 


00 


o 


r^ 


ID 


o 




O 


o 




CM 




CO 


<— H 


00 


LD 


(M 


00 


(M 


1— H 


00 




ID 




CO 


l—i 


O 


ON 




vO 



00 

r- 



CO 



r- 

r-- 



CO 



vO 



CO 



o 

05 

O 

Z 



00 


CM 


ID 


>- 


o 


r- 




CO 




vD 


cr^ 




vO 


ID 




r- 


00 


vO 


ID 


CO 


cn 




O 




CM 


r^ 


CO 


o 


00 


vO 


CO 


00 




1— H 


00 


r- 


ID 




CO 




o 


o 




o 



00 

r- 



co 






CO 



vO 

r- 



W 

Q 

O 



CM 



ID 



vO 



00 






63 



CPU TIME; 56.23 sec 23.50 sec 22.74 sec 
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TEST CASE TWO PROFILE 
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Figure 11 
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RESULTS OF TEST CASE TWO 

















o 


m 


CO 


00 


in 


CO 


CO 


o 


o 


0 




m 




00 


00 


CO 


CO 




00 


CO 


in 


o 


1— 1 


00 




A 




m 


00 


00 




sO 


00 




r— 1 


o 




00 


r- 


in 








ON 


00 


00 


00 


00 


• 


u 


• 


. 


, 


• 


« 


• 




Z 


CO 


00 


CO 


CO 


00 


CO 



















u 




in 


00 




vO 


in 


o 


0 


0 


r— 1 




CO 


f— 1 


vO 


vO 


CO 


yH 


CO 


00 


vO 


00 


in 


00 




A 


CO 




r- 


vO 


CO 




CO 




r— 1 


o 


o 


00 




in 


r- 




ON 




00 


00 


00 


00 


. 


U 


• 


« 


• 


« 


« 


• 




Z 


CO 


CO 


CO 


00 


00 


CO 





O 

S 

Pi 

O 

z 



IT) 

00 

CO 

CO 

I—* 

00 



CT^ 

CO 

O 

CO 



o 

00 

o 

00 

00 



LO 

m 

vO 

00 

00 

00 



00 

00 

00 



h- 

CO 

vO 

LO 

00 

00 



a 

<l> 

CO 

CO 

r-* 

CO 



=3 0 

5 - 


Tf 


CO 


CO 


o 


CO 


00 




•-H 






in 


t^ 




S w 
^ o 


00 


CO 


00 


vO 


in 


00 


CO 




r^ 


o 


CO 


Tf 


f-H 


o 




00 




in 


w g 




o 


00 


00 


00 


00 


« 


• 


• 


• 


• 


• 


z 


00 


00 


CO 


CO 


CO 


CO 



W 

• § 

2 



CO cO in vO 



o 

6 

• r-l 

H 

P 

PU 

O 



CO 

U 

(D 

-M 

0) 

a 

c 

• iH 

CO 

u 

o 

rO 

E 

c 

0) 

> 



oJ 

4j 

c 

O 

N 

O 



65 



TABLE III. 



C. TEST CASE THREE - SYMMETRIC WAVE DUCT 



A profile was constructed corresponding to a symmetric wave 
duct where the sound speed is given by 




a^oJ^CZo-zf 



(115) 



between two depths. The solution to this profile is given by Tolstoy 
and Clay (1966) as 

_ cora.(2w-l)j^^ (116) 

The actual profile dimensions, shown in figure 12, are taken from 
Tolstoy and Clay. The profile was tested at 30.0 and 60.0 hertz. 

The computed eigenvalues agree with the anlytic solutions to one 
part in 10^, with the exception of the higher modes of NORMOl. 

Table 4 shows comparisons between the analytic solution and 
NORMOl, N0RM04 or N0RM05 computed values. 

Of significant interest is the fact that the errors appeared rela- 
tively stable for N0RM04 and N0RM05. The errors for N0RM04 
are consistently 0. 45x10 ^ to 0.69x10 ^meters lowfor 30 Hz, and 
0. 10x10“^ to 0. 12x10“^ meters ^ low for 60 Hz. In contrast, the 
errors for NORMOl increased from 0. 39x10 ^ meters low to 
12. 5x10 ^ meters high for 30 Hz, and from 0. 80x10 ^ to 2. 8x10 ^ 
meters low for 60 Hz. Figure 14 displays a graph of the program 
errors. Figure 13 shows the first three modes. 
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TEST CASE THREE PROFILE 
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TEST CASE THREE MODES 



MODE 1 MODE 2? hOCE 3 




Figure 13 
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TEST CASE THREE ERRORS 



ia.5i X •o*' 




* -K NORM04 and N0RM05 
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.246988417783 0.2469875 - 0.88 0.2469873 - 1.08 0.2469873 - 1.08 
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TABLE 4 




I 




I 




4 



# 





D. TEST CASE FOUR - DEEP SOUND CHANNEL 



Test case four is a typical deep ocean sound channel. The pro- 
file is characterized by a deep and sharp sound channel and by a 
bottom of the same density as the water (Figure 15). With the excep- 

4 

tion of the first mode, the eigenvalues agree within one part in 10 
(Table 5). Note that for the higher frequency, 60 Hz, the agreement 
is somewhat improved. It is interesting that the first mode should 
show a larger disagreement among the programs. This effect may 
be a result of the combination of the sharp sound channel axis and a 
grid spacing of 8 meters. 
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RESULTS OF TEST CASE FOUR 



Frequency: 30 hertz 




DIFF 




DIFF 


MODE 


NORMOl 


N0RM04 


xlO^ 


N0RM05 


xlO^ 


1 


0. 12614128 0. 


12610844 


-32. 84 


0. 12610844 


-32. 84 


2 


567736 


567933 


- 1.97 


567933 


- 1. 79 


3 


537387 


538187 


8. 00 


538187 


8. 00 


4 


514437 


513923 


- 5. 14 


513923 


- 5. 14 


5 


492671 


492224 


- 4.47 


492224 


- 4. 47 


6 


472443 


472294 


- 1.49 


472294 


- 1.49 


7 


453745 


453659 


- 0. 86 


453660 


- 0.85 


8 


436162 


436012 


- 1.50 


436018 


- 1.44 


9 


419601 


419117 


- 4. 84 


419155 


- 4. 46 


10 


404106 


400638 


-34. 68 


402169 


-19. 37 


Frequency: 60 hertz 










1 


0.25255354 0. 


25250819 


-45. 35 


0. 25250819 


-45.35 


2 


196659 


196827 


1. 68 


196827 


1. 68 


3 


154658 . 


154783 


1. 25 


154783 


1.25 


4 


118755 


118986 


2. 31 


1 18986 


2.31 


5 


089177 


089755 


5. 78- 


089755 


5.78 


6 


063654 


063571 


- 0. 83 


063571 


- 0.83 


7 


039561 


039343 


- 2. 18 


039343 


- 2. 18 


8 


016594 


016618 


0. 24 


016618 


0. 24 


9 


0.24994761 0. 


24994890 


1. 29 


0. 24994890 


1.29 


10 


973743 


974255 


5. 12 


974255 


5. 12 



TABLE V 
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E. TEST CASE FIVE - DOUBLE CHANNEL 



Test case five is a double sound channel of nearly symmetric 
dimensions (Figure 16). This profile was used to test the multi-duct 
properties of the N0RM04 and N0RM05 programs. Note that for 
30 Hz thediffer ences between NORMOl and N0RM04 (or N0RM05) 
increase by almost an order of magnitude for modes 5 and 7 (Table 6). 
For these modes the barrier connection formula (equation 114) is 
employed and the ducts are considered ’’connected. The mode pro- 
files (Figure 17) for these modes also do not completely correspond. 
Modes 5 and 7 appear to be ’’upside-down” in N0RM04, with the 
dominant amplitude waves in the lower channel. This is due to the 
fact that N0RM04 and N0RM05 iterate the vertical wave number 
from the surface downward. Because of this, the phase at the upper 
turning point of the lower channel represents the phase angle for the 
third and fourth mode in the lower channel. Equation (114) remains 
correct, although here it has been incorrectly applied. 
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RESULTS OF TEST CASE FIVE 









DIFF 




DIFF 

/ 


MODE 


; NORMOl N0RM04 


xlO^ 


N0RM05 


xlO® 


Frequency: 30 hertz 








1 


0. 12608464 


0. 12604450 


-40. 14 


0. 12604450 


-40. 14 


2 


08403 


04450 


-39. 53 


04450 


-39. 53 


3 


0. 12553486 


0. 12554660 


11. 76 


0. 12554660 


11. 76 


4 


52366 


53937 


15. 71 


53937 


15.71 


5 


16376 


36343 


199. 67 


36343 


199. 67 


6 


11170 


09716 


-14. 54 


09716 


-14. 54 


7 


0. 12485241 


0. 12490245 


50. 04 


0. 12490245 


50. 04 


8 


75394 


74793 


- 6. 01 


74793 


- 6. 01 


9 


59913 


59806 


- 1.07 


59806 


- 1. 07 


0 


42282 


42565 


2. 83 


42568 


2.86 


Frequency; 60 hertz 








1 


0.25248206 


0.25242998 


-52. 08 


0. 25242998 


-52/08 


2 


48204 


42998 


-52. 06 


42998 


-52. 06 


3 


0.25179326 


0. 25180252 


9. 26 


0. 25180252 


9.26 


4 


79260 


80252 


9. 92 


80252 


9.92 


5 


31541 


31134 


- 4. 07 


31134 


- 4. 07 


6 


30922 


31100 


1. 78 


31100 


1. 78 


7 


0.25087788 


0.25088532 


7. 44 


0. 25088532 


7.44 


8 


84379 


84790 


4. 11 


84790 


4. 11 


9 


50385 


49918 


- 4. 67 


49918 


- 4. 67 


.0 


41349 


40770 


- 5. 79 


40770 


- 5.79 






TABLE VI 
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DOUBLE CHANNEL MODE COMPARISONS 
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% 







VI. CONCLUSIONS 



A. GENERAL 

As seen by table 7, the WKB method has performed consistently 
faster than the iterative technique. The WKB programs (N0RM04 
and N0RM05) found solutions with about one -half the Central Pro- 
cessing Unit (CPU) time required by the finite difference program 
(NORMOl). In addition, one should consider the method used for in- 
tegration of the vertical wave number in the WKB programs. The 
trapezoid rule was used which required the evaluation of a square 
root at each grid point. Although simple, this is a time consuming 
task, while it is not necessarily very accurate. If a faster integra- 
ting procedure were substituted for the trapezoid rule, the time re- 
quired for each iteration could be reduced considerably. Thus, in 
addition to showing an improvement in computation time, the WKB 
method also offers the potential of greater computational speed 
through the use of sophisticated integrators. 

The WKB and finite difference methods appear to have similar 
accuracy in comparison with other programs. The accuracy attain- 
able is adequate to calculate useful modes for a propagation loss pro- 
gram. However, the errors involved in the differences between 
successive eigenvalues will not permit accurate coherent phasing of 
the modes beyond the first or second convergence zone . Coherent 
phasing at extremely long ranges is, in fact, attainable only with an 
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CPU TIMES FOR VARIOUS RUNS 





NORMOl 


N0RM04 


N0RM05 


Test Case One 


56. 23 sec 


23. 50 sec 


22. 74 se 


Test Case Two 


27.43 


14. 73 


14. 92 


Test Case Three 


106. 35 


38.42 


39.47 


Test Case Four 


105. 67 


40. 23 


40.43 


Test Case Five 


114. 07 


44.53 


47. 09 


TOTAL 


409. 75 sec 


161.41 sec 


164. 65 se 



TABLE VII 



analytic profile and solution — an ideal mathematical model. In 
fact, when one considers the accuracy of the available sound velocity 
data, it is difficult to imagine any model attempting to yield accurate 
phasing effects at any large distance. Since the position of a target 
is unknown, in practice it is more important that the amplitude and 
nature of phasing effects be know, than for the sharp interference 
peaks to be accurately predicted. Thus, accuracies of one part in 

5 

10 should enable a model both to take account of the phasing effects 
in the first convergence zone, and to give a general description of the 
effects of phase interference at longer ranges. 

Although the WKB programs do not correctly interpret the barrier 
effects, the problem appears to lie in the approach taken in the pro- 
gramming, rather than in the mathematical and physical relationships. 
Within any duct, the effects of any upper and lower ducts must be 
treated as a turning point phase effect similar to the top and bottom 
boundary conditions. The modes 5 and 7 of the N0RM04 30 Hz double 
channel case are unwitting examples of this. In N0RM04 and 
N0RM05 the ''connected** ducts are considered as a single wave sys- 
tem. In fact, with any sized barrier, no matter how small, one must 
consider the wave system within each duct separately. In this case 
the turning point phase angle between the duct and barrier is analo- 
gous to the impedence due to the other duct felt by the wave system in 
its duct. 
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B. FUTURE PROGRAM 



It appears that a faster WKB method normal mode program can 
be written which will properly interpret a barrier. In order to pro- 
perly consider the multiple channel cases, a mapping procedure 
should be used before actual mode searches commence. The sound 
profile would be divided into "levels” (Figure 18 ), found by selecting 
all the velocity maxima and minima. Each "level" represents a 
range of the possible eigenvalues, Within each "level" there is 

a single combination of ducts, barriers and boundaries. At the start 
of computations the program would determine the limits of each level. 
Then, at the beginning of the frequency loop, the program would de- 
termine the maximum and minimum number of modes for each duct 
within each level. Thus, equipped with a map of the profile charac- 
teristics, the program would compute eigenvalues for each duct and 
level in sequence (Figure 19 ). After all the eigenvalues are computed 
and sorted (some may be out of numerical order), the eigenfunctions 
can be computed quickly using the finite difference procedure of 
NORMOl. Combined with a faster integrating method, such a pro- 
gram should yield fast, accurate normal modes for an arbitrary 
sound speed profile. 
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depth 



SOUND VELOCITY PROFILE 
WITH 
LEVELS 




Figure 18 



83 



DEPTH 



SEQUENCE OF MODE CAECUEATIONS 




b Maximum Modes Region Z 
c Maximum Modes Region 3 
d Maximum Modes Region 4 



Figure 19 
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APPENDIX A 



NORMOl FINITE DIFFERENCE SCHEME 



In order to iterate equation (18) through the depth grid, a finite 
difference scheme is used which computes values of the eigenfunction, 
Z(z) and its derivative, Z'(z) at each grid point. This scheme was 
originally developed by the Arthur D. Little Company [Report No. C- 
70673, 31 January 1969]. 

Z(z) is expanded into a Taylor series, 

=Zr^) * hZm * h!Z''z) - , (91) 

2! 3! 

where h is the depth grid spacing. From equation (18) we have 

, = - [E"\/(h)] Z< 2 ) = - F(z)Z(z) , (92) 

where F(z) represents the square of the vertical v/ave number. Using 
a forward difference for Z"'(z), 

2% = (93) 

h 

Substituting equation (92), we have 



H^Tz) = i[F(z)2u.) ~ F(z*h)Z(z*h)J , 

h 

Combining equations (91), (92) and (94), 

2? 3*[ h 



(94) 



(95) 



85 



we have on clearing, 



-- ZfE)]\-h"F(H)j ^ hZ'. 



(96) 



Placing the above result in a vector subscript notation we have the 
formula used to compute the succeeding eigenfunction, 



Zui = Zc(i-^Vi)v hz"c 



(97) 



We must now compute the value of the derivative of the eigenfunction 
Z*(z) at the next grid point. In order to do this let us express the 
eigenfunction at our original grid point as a backwards Taylor series 
of the new grid point, z+h, 



2? V. 



(98) 



Substituting equations (92) and (93), we have 

Z (r) - h) - b.^ F(£.-W')Z( zi>k) - ^ (99) 

Rearranging, we have, 

s P(lfh)^ “ Z(E>{l ^ j 



( 100 ) 



or in vector notation, 

H'l.i = 1 [Zui (1- Fc.t) - Hi. U * g FO] 



( 101 ) 
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APPENDIX B 



N0RM04 CONNECTION FORMULAE 



N0RM04 requires the evaluation of the phase at the upper and 
lower turning points in addition to the phase effect of a barrier. For 
the purposes of this discussion we will align the origin of the depth 
axis with the turning point, designating the positive direction towards 
the ’*real’* region. 

Consider the upper turning point case (Figure 5). The free sur- 
face boundary condition requires that the "imaginary’’ region solution 
disappear at the surface. From equation (59) this requires that 






( 102 ) 



Thus, we can define C and D (except for a constant of multiplication); 



C-- -e 



-St(-H) 



D = e 



Sz(-H) 



( 103 ) 



Applying equation (71), we have 



^ _ gSz(-H) 

"b ” g-Sa(- n ) g St(- h) 



( 104 ) 



or 



~ Arc ban |^Tan h 



( 105 ) 
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Now, consider the bottom boundary condition and lower turning 
point (Figure 5). The bottom boundary condition (equation 26) re - 
quires 



^Kb = Ks 



( 106 ) 



where is defined as the imaginary vertical wave number evaluated 
s 

at the water side of the bottom interface, 



C(-H) 



(107) 



Thus, the imaginary solution can be written by defining C and D such 
that 



C= (PbKs " ^oKb) e 
D= (fbKs- 



Now, by applying the turning point connection formulae (equation 71) 
we have 



g. (ebKs* ^ 

‘^'(ebKs-eoKB)e^*^-”) - (^bKs ^ ’ (109) 



or 



©L = Arc tan 



D^ -1 






( 110 ) 
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where is defined as 



Da- ebKs - eoKs 
CbKs CoKr 



( 111 ) 



Now, consider the situation with an upper and lower channel 
separated by a barrier (Figure 8). Let the z origin be at the lower 
connection point, z^, with the positive direction downward. From 
equations (59), (69) and (71) we have for a connection at the upper 
point, Zj, 



ai _ De~^ Ce*- 
bi De-'- - Ce*- ’ 



( 112 ) 



where L. is as defined in equation (86). Letting Aj equal the numer- 
ator and the denominator, we have, except for a constant, 



C=|(a^-bOe-‘- 

D=|(ai , 



(113) 



Applying these expressions to equation (71), we have 



o-z (ai + bi)e‘- 



(114) 



(cLi + bi)e‘- - (ai- bi.)e 



-u 



or 



fl, = A.,K.^r (tat,ei-l)e-.aan9i-l)e-^~ | . (115) 

[_ ( to.n 9 x ♦ 1 ) e*" - ( tan 01 - 1 ) e ‘"J 
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APPENDIX C 



N0RM05 CONNECTION FORMULAE 



For N0RM05 the development of the upper and lower ’’imaginary” 
region particular solutions by the application of the free surface and 
bottom boundary conditions is identical to that for N0RM04. The 
differences in the two programs lie in the application of the Airy 
phase connection formulae vice the asymptotic connection formulae. 

To obtain the upper turning point phase, we apply equation (103) to 
equation (77) 



a . ±(Z 



g^-SzC-H) J 



(116) 



b = 



/I 






(117) 



or 



Gu = Arctan 



rD±-ii 

LDi.1_ 



where is defined as 



Di = 



(118) 



(119) 



Note that as H approaches zero (the turning point approaches the sur- 
face) the value of 9^ in equation (119) approaches arctan (1/3) vice the 
phase we required for the free surface boundary condition, (0^ = 0 ). 
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To derive the lower turning point phase, we apply the particular 
lower ^hmaginary region solution (equation 108) to the Airy phase 
connection formulae (equation 77), 



a - i r2(et,Ks- . (ebKs* 



b' =|r [2(^1. Ks - ec Kb) 



( 120 ) 



( 121 ) 



Dividing, we have 



0^ = A«ta. [f^] , 



( 122 ) 



where is defined as in equation (HI). Note that this equation also 
does not approach the bottom reflected boundary condition (equation 
84) as H approaches zero (the lower turning point approaches the 
bottom). 

Lauvstad (1971) derived an expression for the phase coefficients 
of a wave entering a barrier region in terms of the phase of the wave 
leaving the opposite side of the barrier: 



\(o.z ^ 

L 



(CL^- 





(OLz. - 




( 123 ) 



By rearranging, we obtain 
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(124) 



be - (at-bi)e'^J 

which agrees with equation (114), the N0RM04 result, 

Gl = Archan Rai ^ U)e‘- ^ (ai- bQe'^ 

I (at^-bt)e‘- - (ai-bi)e'^ 



(114) 



APPENDIX D 



INPUT DATA 

The following format is used for the input data to all three programs, 
NORMOl, N0RM04 and N0RM05. 



VARIABLE 


FORMAT 


MEANING 


RANGE 


CARD ONE 








HED (20) 


20A4 


HEADING. TITLE CARD 




CARD TWO 








NUMV 


no 


Number of Depth/Sound Speed 
Pairs. Negative if conver- 
sion from feet to meters de- ’ 
sired. 


2 - 50 


NMOD 


no 


Number of Modes Desired 


1 - 100 


lEX 


no 


Iteration limit. 


1 - 14 


N 


no 


Number of grid points desired 


o o 

0 o 

1 

' 

o 

o o 
o o 


lOUT 


no 


Type of output desired 


see below 


IDEV* 


12 


Output device for cardfile 




CARD THREE 








H 


Flo. 3 


Bottom Depth 




CB 


FIO. 3 


Bottom Sound Speed 




DRO 


FIO. 3 


Water density 




DRB 


FIO. 3 


Bottom density 




HP LOT 


FIO. 3 


Maximum depth to be plotted 


HPLOT H 



*N0RM04 and N0RM05 only. 
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CARD FOUR 



NFREQ 



no 



Number of frequencies 
desired 



CARD FIVE (SIX) 



FREQ (NFREQ) 8F10.3 Frequencies desired 

CARDS SIX TO 5+NUMV (SEVEN to 6+NUMV) 



DEP(I) 
CC (I) 



FIO. 3 
FIO. 3 



Depth 

Sound Speed 



1 - 20 



DEP(I)> 

DEP(I-I) 



For NORMOI, IOUT = I gives a punched deck output of the sound 
profile and mode parameters and profile. For N0RM04 and 
N0RM05, the lOUT parameter gives the following types of output; 



lOUT PRINTER GRAPH PUNCHED CARDS DISK OUTPUT 



0 


YES 


no 


no 


1 


YES 


YES 


no 


2 


YES 


no 


YES 


3 


no 


no 


YES 


4 


no 


YES 


no 


5 


no 


no 


no 
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NORMOl FLOWCHART 
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FREQUENCY LOOP 





SZ 

USE METHOD OF HALVES 
TO GET UPPER AND LOW- 
ER BOUND OF EIGENVALUE 
WITH M-1 AND M MODE 
CROSSINGS 



C 



:sz 

ITERATION LOOP BEGINS 



Ph 

o 

3 

w 

o 



) 



^ 


7 


COMPUTE r 
BY y 


■lEXT GUESS 
ALVING 








INCREMENTAL 
CHANGE < DESIRED 
ACCURACY? 




yes 



no 



CALL ITERAT 






HERAT lOi; 
COUNTER 
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I 











FREQUENCY LOOP 



15 12 4 



P4 

o 

S 



w 

Q 

S 







10 



CHECK FOR DEGENERATE 
SOLUTION AND FIXUP 
IF REQUIRED 



.. X 

^MODE HAS BEEN FOUN^ 



I 



NORMALIZE MODE TO ONE 






zi 





FREQUECNY LOOP 
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N0RM04 AND N0RM05 FLOWCHART 




99 



FREQUENCY LOOP 




100 



FREQUENCY LOOP 




101 



SEARCH SUBROUTINE FLOWCHART 




102 




103 




20 




104 



i 




I 



I 

I 
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WKB SUBROUTINE FLOWCHART 




106 
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O' 




CCCCCCCCCCCCCCCCCCCCCC.CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



c c 

C NCRMOl : A NORMAL MODE EIGENVALLE C 

C AND EIGENFUNCTION SOLVER C 

c c 

C LT KIRK E. EVANS, USN C 

c c 



ccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

IMPLICIT REAL*8(0) 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLOCKS >>>> 
C 

COMMON /DBLL/ DV Z ( 50 1 ) , CH , DH2 , DC&S G ,ORO RB , CKAP , DUMAK 
COMMON /DBLE2/ D , OCB , DC M I N , DW , CC L , DC 1 S G 
COMMON /D6LE3/ DUZ ( 50 1 ) ,OZZ ( 50 I ) 

COMMON /DBLE4/ D S , CSS CC , DP V EL , D G V E L 
COMMON /PARAM/ N,NP1,IEX 
COMMON /SOUND/ Z(50),C(50J 
COMMON /HEAD/ H.ED(ZO) 

COMMON /SINGl/ ZM , C8 , CM IN , CMAX , PLTMAX ,UM , FG 
COMMON /INPUT/ NUMV,NMOO 
COMMON /DMAP/ OKI (LOO) 

COMMON /DKMAP/ DKM IN ,OKMAX , DKL , DKU 
COMMON /DSGUNO/ 0CC(501) 

COMMON /RHC/ RO,RB 

COMMON /BOTTO/ DEC AY , ZB ( 20 ) ,0B ( 20 ) 

COMMON /FPEG/ FQVI20) 

COMMON /ACC/ DEPS,OIFFI 
C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DATA >>>»»>>»» 
C 

DATA DP2/6. 28318530600/ 

DATA NPRINT/6/ 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS >>> 
C 

c 

WRITE (NPRINT,16) 

DC 1 1=1,100 

1 DKK I) =-l ,CDO 
C 

CALL INPUTl (IOUT,NFREG) 

CALL INTRPL 
DCMIN2=DCMIN=?'DCMIN 
CALL CZFLQT 
IF (lOUT.NE.l) GO TO 2 
CALL PRCFIL ( CZZ , DCC , NP 1 ) 

2 CONTINUE 
NMODP=NMOD 

C 

C . FREGUENCY LCCP BEGINS 

C 

C 

CC 15 IF=1,NFREQ 
FC=FQV ( if; 

DW=DP2"FG 

DWZ=DW-<'DW 
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i 





oooo ooooo o o o ooo o o o o ooo oooo ooc^ oooo 



DKAP=DW/DCMIN 



-COMPUTE THE V(Z) FUNCTION 



DC 3 

CCZ=DCC ( J J 
DCZ2=CCZ*DCZ 

2 DVZ(J)=DW2=^=( (CCZ-DCMIN)*(OCZ + DCMIN) ) / ( DC Z2 =^CCM IN2 ) 



WRITE (NPRINT,18) FQ 

NPCD=NMCDP 

CALL MAXMIN (£13) 

mode loop begins 



DC 12 M=1,NM0D 
IT=1 

CALL HALF (M,G9) 

ITERATICN LOOP EEGINS 

A DKN=( DKL+DKU)*0.5D0 

DIKC=DABS(DKU-DKL) «C.5DC/DKN 
IF (DINC-OIFFI) 10,5,5 

5 CALL HERAT tOKN.MS) 

IT=IT+1 

DEPSIL=DUMAX*DEPS 

IF (DABS{CUZ(NP1) )-DEPSIL) 11,6,6 

6 IF (MS-M) 7,8,8 

7 DKU=DKW 
GG TO 4 

8 CKL=DKN 
GC TO 4 

MODE HAS BEEN FCLND 

9 DKN=( DKL+DKU)*0.5DC 
1C CALL FIXUP (DKN) 

11 DKMAX=DKN 
DK1(M)=DKN 

CALL DNCRMl ( DUZ , DUMAX , NPl ) 

CALL BOTTOM (DKN,M) 

CALL MODPLT (M, I ) 

WRITE (NPRINT,17) IT 
CALL INTEGR (DKN,M) 

IF (lOUT.NE.l) GO TO 12 
CALL OUTPUT ( DUZ , DKN , M ,NP1 ) 

12 CONTINUE 



PCDE LOOP ENDS 



12 DC 14 11=1,100 

14 DKK II ) =-l .000 

15 CC^TINUE 



FREQUENCY LCCP ENDS 
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ooo 



STOP 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 



It FORMAT Cl K.E. EVANS NORMAL MODE EIGENVALUE', 

=?--'SCLVER, VERSION 1, REVISION 6, DATED 13 AUGUST 73') 

17 FORMAT CO NUMBER OF ITERATIONS 15) 

18 FORMAT CO FREQUENCY : F6.1, • HZ.') 

END 
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SL6R0UTINE INPUTl (lOUT^NFREQ) 

INPLICIT REAL=^8(D) 

c =*< 

C * THIS SUBROUTINE INPUTS THE SOUND AND PIN DATA, 

C * DISPLAYS IT, THEN SETS SOME CONSTANTS FCR 

C * LATER USE. 

C * 

c 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COPPCN BLOCKS »»»>»» 
C 

CC^MON /DBLl/ DVZI501) ,DH,DH2,OCBSC,DRORE,CKAP,DUF'AX 
COMMON /DBLE2/ D W2 , DC B , DC MI N , D W , DC 1 , DC I SC 
COMMON /HEAD/ HEC(2CJ 

COMMON /SINGl/ Z M , 0 6 , C MI N , C MAX , P LT MAX , UM , F C 

COMMON /SOUND/ Z(50),C{50J 

COMMON /INPUT/ NUMV,NMOD 

COMMON /PARAM/ N,NP1,IEX 

COMMON /DHhH/ DH2S IX, DH2D3 , DHHLF 

COMMON /RHC/ RO,R6 

COMMON /FREQ/ FQV(20) 

COMMON /ACC/ DEPS,DIFFI 
C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DATA >>>>>>»»>>>» 
C 

DATA NREAD ,NPRINT/5,6/ 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS >>>» 

C 

C 

CCNVRT=1.0E0 

C 

Cl READ IN TFE RUN PARAMETERS 

C 

READ (NREAD, 4) HED 
V.RITE (NPRINT,5) HED 
C 

READ (NREAD, 6) NUM V , NMOD , I EX, N , I CUT 
N=IABS( NI 

IF (N.GT.500) N=50C 
IF (N.LT.5C) N=LOO 
NMCD=I ABS (NMOD) 

IF (NMOD. GT. LOO) NM0D=100 
IF (NMOD. EC. 0 1 NMGC=10 
IF (N.EC.OJ N=500 
IEX=IABS( lEX) 

IF (IEX.GT.L4) IEX=I4 
IF ( lEX.EQ.O) IEX=7 
NFl=N-H 

READ (NREAD, 7) ZM , CB , RO , R B , PLTMAX 
ZM=ABS( ZM) 

PLTMAX=ABS(PLTMAX) 

READ (NREAD. 6) NFREC 
NFPEQ=IA6S(NFREQ) 

IF (NFREC. GT. 20) NFPE0=20 

READ (NREAD, 7) ( FG V ( I ) , I = 1 , NFRE C ) 

FC=FCV( 1) 

C 

C— . ^-CONVERT FPCM FEET TC METERS 

C IF REQUIRED 

C 

IF (NUMV) 1,2,2 

1 CCNVRT=0.3048 
NUM V=-NUMV 
ZM=ZM*COMVRT 
Ce=CB«CCNVRT 
PLTMAX=PLTMAX*CONVRT 

C 

2 WRITE (NPRINT,8) 

WRITE (NPRINT,9) ( FQV ( I ) , I = 1 , NF REC ) 
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oooooo o oooo oo oooo 



WRITE (NPRINT,10) NPOD.IEX 
WRITE (NPRINTrll) ZP,PLTMAX 
WRITE (KPRINT,12) R0,R6 
WRITE (KPRINT,13) CB,N 



RE/yQ DISPLAY SCUNO DATA 

WRITE (NPRINT,14) h£D 



DC 3 I=1,NUMV 

READ (NREAD,7) Z(I),CI 

CI=CI’i'CCNVRT 

Z( I ) = Z( li’f^CONVRT 

C ( I ) = CI 

3 WRITE (NPRINT,15) Z(I) ,CI 



COMPUTE SCRE CONSTANTS 

Dh = CBLE(ZM)/DFLOAT (i\) 

D2R=D6lE('ZM) 

CRCRb=DELE(RO/RB) 

CR2=DH^‘CH 

CCe=D6LE(CEI 

ccesQ=cce=!=cc8 

DE2SIX=DH2^0. I 66666666 6666667D0 
DR2D3 = Dh2*«0.33333 33333333 333D0 
DhRLF=Dh^0.5D0 
DEPS= 1. CDpi'-^I-IEX) 

CIFFI=D£PS^DEPS 

CIFFI = DRAXUDIFF I ,1.0D-16) 

RETURN 



<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 



A FORMAT 

5 FCRMAT 

6 FCRMAT 

7 FCRMAT 

8 FCRMAT 

9 FORMAT 
1C FORMAT 

1 

11 FCRMAT 
1 

12 FCRMAT 

* •. 

1 

13 FORMAT 

* F8 

I 

1^ FCRMAT 
1 

15 FCRMAT 
END 



(20A4) 

{• 1' ,20A4) 

(5110) 

(SF10.3J 

(•O', T34, 'RUN PARAMETERS') 

( • O' ,T30, ' FRECUENCIES =', (T42,F6.l,' HZ',/)) 
(•0',T24, 'MODES REC.UESTEC :', 14, //, 

T24, • ITERATION LIMIT = IO“=s^( -',12,')') 
('C',T27, 'BOTTOM DEPTH :',Fe.l,' METERS',//, 

T22, 'MAX DEPTH PLOTTED :',F6.1,' METERS') 
C0',T26, 'CCEAN DENSITY :', F7.4, 

GM/CM^>i^3', //, 

T25, 'BOTTOM DENSITY :',F7.4,‘ GM/CM=»=*3') 
C0',T18, 'BCTTCM SOUND VELOCITY 
1,:' METERS/SEC',//, 

T17, 'NUMBER OF GRID POINTS :', 15) 

Cl', 20A4, //, T33, 'VELOCITY PRCF I LE • ,// , T22 , 
•DEPTH, METERS' ,T44, 'SOUND VELOC I TY ,M/ScC ' , / ) 
('O' , T25, F6.1, T57, F6.1 ) 
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SieROUTINE IKTRPL 
INPLICIT REAL*8(D) 






SOUND 



THIS SUBROUTINE INTERPOLATES THE DEPTH /iND 
VELOCITY DATA GIVEN, FOR THE GRID VALUES. 

IN ADDITION, IT CCiiVPUTES SOME CCNSTANTS FCR 
LATER USE, AND EXTRAPOLATES THE SOUND VELOCITY 
TO THE BOTTOM IF RECUIRED. 



CCPMCN /DBLl/ DVZ(50U ,CH,DH2,DCeSC,CR0RE,CKAP2,CUP.AX 
COMMON /DELE2/ DW2 , DC B , DC M I N , Dt. , DC I , DC 1 S C 
COMMON /DDLE3/ OUZ ( 501 ) , D2Z ( 50 1 ) 

COMMON /PARAM/ N,NP1,IEX 

COMMON /SINGl/ ZM , CB , CM I N , C MAX , PLTMAX , UM , FG 
CCMMON /USCUNC/ DCC(501j 
COMMON /SOUND/ Z(5C),C(50) 

CCMMON /INPUT/ NUMV,NMOC 



CMAX=0.0 

CMIN=CB 

EXTRAPOLATE THE PROFILE TO THE BOTTOM *** 
IF (Z(NUMVJ-ZMJ 1,2,2 

1 NLMV=NUMV+1 

2 C (NUMV)=C(NUMV-l)-+(ZM-Z(NUMV-l) )*C. 17E- 1 
Z (NUMV) =ZM 

3 NUMVPl=NUMV+l 
NLMVD2=NUMV/2 



<<<<<<<<<<<<<<<<<<<<<<<<<<<< SWAPPING AND TESTING LOOP 

DC 4 I=l,NUMVI>2 
NLF=NUMVPL-^I 
CL=C( I ) 

CU=C(NUP) 

ZL = Zl I ) 

ZL=Z(NUP) 

CMAX=AMAXUCMAX,CL,CU ) 

CMIN=AMINI (CMIN, CL,CU) 

C(I)=CU 
C(NUP)=CL 
Z( I) = ZM-ZU 
Z{NUP)=ZM-ZL 
^ CONTINUE 

<<<<<<<<<<<<<<<<<<<<<<<<<<<< SWAPPING AND TESTING LOOP ENOS 

IF (NUMVD2^2.EQ.NUNV) GO TO 5 
CNV021=C ( NUMVD2+1 } 

CMAX=AMAXUCMAX,CNVC2l ) 

CM IN=AM IiN UCM IN, CN V02 1) 

Z(NUMVD2+l )=ZM-Z(NUMVD2+1 ) 

SET UP FQR TTIE FINE GRID COMPUTATION ’K’** 

E DCMIN=D8LE (CMIN) 

CCMIN2 = CCMIN>--0CMITJ 
DC1=DBLE(C(1) ) 

CC1SQ=CC I’f'DCl 
C2M=0BLE(ZM) 

I = l 

CCI=DC 1 
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oooo ooo ooo ooo oooo 



DCIP1=CEL6(C(2) ) 
D2I = 0BLE( Z(U J 
C2IPl=CeLE(Z(2) J 
02DIFF=DZI Fl-DZI 
CFLN=DFLOAT(N ) 
C2F=C.0C0 



<<<<<<<<<<<<<<<<<<<<<<<<<<<< FINE GRID LOOP BEGINS 



DC 8 J=l,NPl 
D2Z( J)=CZP 

IF (DZP-DZIPL) 7,7,6 

RESET THE I PARAMETERS 



6 I=I+L 
D2I=0ZI FI 
CCI=DCIP1 

CCIPi = CELE(C( I+Ii ) 

CZIP1 = DEL5(Z(1 + U ) 

C2CIFF=0ZIP1-0ZI 

COMPUTE THE VALUES FOR C(2) *** 

7 DCZ=(DCI^(OZIPL-DZP )+DCIPl*(DZP-CZI ) )/DZDIFF 
CCC( J) =DCZ 

COMPUTE THE NEXT DEPTH 

C2P = CZM*DFLnAT(J JVDFLN 

8 CONTINUE 



<<<<<<<<<<<<<<<<<<<<<<<<<<<< FINE GRID LOOP ENDS 



RETURN 

ENC 



IL4 



SIBHOUTINE MAXMIN (♦) 
IMPLICIT REAL^8(0) 



C 

c * 

C * THIS SUBROUTINE CHECKS THE NAXIHUM AND MINIMUM 

C « VALUES OF THE HCRIZGNTAL WAVE NUMBER; THEN 

C * DETERMINES THE NUMBER OF MODES REPRESENTED 

C * BY EACH. IF NC MODES ARE PRESENT FOR THE 

C « MINIMUM WAVE NUMBER THE FRECLENCY IS BELCW 

C ^ CUT-CFF. IF THE NUMBER OF MODES REQUESTED IS 

C * ABOVE THAT REPRESENTED BY THE MINIMUM ^AVE 

C * NUMBERf THE REQUEST IS MODIFIED ACCORDINGLY. 

C * 

c 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLOCKS >>>>>> 
C 

c 



COMMON /DBLl/ OVZ ( 501 ), DH , 0H2 , OCBSG , DROP B , DKAP , DUMAX 

COMMON /DBLE2/ DW 2 , DC B , CCM I N , DW , DC i ,DC I SQ 

COMMON /DBLE3/ DU Z ( 501 ) , DZZ ( 50 I i 

COMMON /PARAM/ N,NP1,1EX 

COMMON /HEAD/ HED(20) 

COMMON /SINGl/ ZM , CB ,CM IN , C MAX , PLT MAX ,UM , FG 
COMMON /OMAP/ DKl(lOO) 

COMMON /OKMAP/ DKM I N , DKMAX , DKL , OKU 
COMMON /INPUT/ NUMV,NMuD 
COMMON /SOUND/ Z(50),C(50) 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DATA >>>>»>>»>»» 
C 

DATA NREAD ,NPRINT/5 ,6/ 

DATA MAXI,MINI/'MAXI* , 'MINI'/ 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM EEGINS »>» 

C 

C 

C COMPUTE THE MAXIMUM AND MINIMUM DK'S 

C 

CKMIN=OW/DCB 

DKMAX=DW/DCMIN 

C 

C *** CHECK OUT THE OKMIN FOR MAX NO. OF MODES 

C 

CALL ITERAT (DKMIN,KSI 
IF (MS) 1,4,1 
1 MODMAX=MS 
DK1(MS)«DKMIN 
C 

MCCMIN=0 
MS = 0 



C 

C *** PRINT THE RESULTS *** 

C 

WRITE (NPRINT,5) HED 
WRITE (NPRINT,6) 

WRITE (NPRINT,?) MAX I , DKM AX , MOD MI N 
WRITE (NPRINT,?) M I N I , DKM I N , MOD M A X 
C 

C *’!'* CHECK THE NUMBER CF MODES REQUESTED 

C 

IF (MCDMAX-NMOD) 2,3,3 
2 NMCD=MOOMAX 

WRITE (NPRINT, 3) 

WRITE (NPRIfNT,9) 

WRITE (NPRINT, 10) 

WRITE (NPRINT, 9) 

WRITE (NPRINT, 3) 

WRITE (NPRINT, 11) NM.OD 
C 
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2 RETURN 
C 

C *=»>? BELOW CUT-OFF FREQUENCY 

C 

^ WRITE {NPRINT,8) 

WRITE (NPRINT,9> 

WRITE (NPRINTtlO) 

WRITE (NPRINT,9) 

WRITE (NPRINT,8) 

WRITE (NFRINT,12) 

RETURN 1 




FORMAT ST/^TEMENTS » 



IC 

11 

12 



5 FORMAT 

6 FORMAT 
1 

7 FORMAT 
1 

8 FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 



1 

2 



Cl', 20A4) 

CO', ///, T25, 'BOUNDS 
' WAVE NUMBER' J 
CO', T22 , A4 , ' MUM K : ' 
'MODES :', lA) 

IT33, 13 ('*')) 
(T33,1H-,11X,1H«) 

(T33,'=? CAUTION *') 
CO* 

( ' 



ON HORIZONTAL' 
, F12.8, 3X, 



T26, 'MCOES REGUESTEC RESET TC 
. _ , T25, 'FREQUENCY RECUESTED I^ 

BELOW CUT-OFF', //, T20, “ 

NEW FREQUENCY SELECTED.') 



0 



I A) 



'JOB TERMINATED OR • 



END 
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SLEROUTINE HALF ( M , ♦ ) 
I^FLICIT REAL’^-8(D) 



* THIS SUBRCUTINE SEARCHES FOR AN UPPER VALUE 
» OF CK WITH M-1 MODE CROSSINGS AND A LOWER 

=!^ VALUE CF DK WITH M NODE CROSSINGS. AT THE 

=«^ SAME TIME IT FILLS IN A MAP CF DK VS. 

MODE CROSSINGS FOR LATER USE. 

* 

S{: ^ ^ * :{t y: * :^c A * !j! * * :(« :4t * * * * 5^ :{: 5{; jf: 4: 5je ^ 2ic ^ jf: ;Ji jjc ^ :<c # 3{; !(« + :(t ;{[ * 



<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<CCMMCN BLOCKS >>>>>>>>>>>>> 
COMMON /DBLE2/ DW2 t CC B , CCM I N , 0 W , DC I ,OC I SC 
COMMON /PARAM/ N,NP1,IEX 
COMMON /DMAP/ OKU ICO) 

COMMON /DKMAP/ DKM I N , DKMAX , CKL , DKU 
COMMON /ACC/ DEPS,CIFFI 



<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< program begins >>>>>>>>>> 



LOW 

LUF 



I 

1 



**=5“ CHECK THE MAP FOR ALREADY COMPUTED GOESES **« 



DC 2 I=M, 100 
IF (DKUI) ) 1,1,3 

1 CONTINUE 

2 CONTINUE 



DKL=DKMIN 
GO TO 4 

2 CKL=DK1( I) 
A CKU=DKMAX 



FIND C CHECK THE NEW MID DK VALUE *** 

5 CKN = ( DKL+DKU)=!'0.5D0 

*** CHECK THE DIFFERENCE 

IF (( CKU-DKN)/DKN-DIFF I ) 13,6,6 

6 CALL ITERAT (OKN,MS) 

IF (MS-M) 7,8,9 

*** IF WE GET TWO END POINTS WE ARE FINISHED 
*** IF NOT, CONTINUE HALVING 

7 CKU=DKN 

IF (LUP.GT.OJ RETURN 
LCW = 1 
GC TO 5 



8 CKL=DKN 

IF (LOW.GT.O) RETURN 

LUF=1 

GC TO 5 

FILL IN THE MAP IF NEEDED 

9 IF (MS-ICC) 10,10,6 

1C IF (CKN-DKKMS) ) 12,12,11 
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11 CK1(MS)=0KK 

12 CKL=OKN 
GC TO 5 

* + ^ IF WE GET HERE WE HAVE COME TC WITHIN =(■>*« 
*** THE LIMITS OF DESIRED ACCURACY IN DK 

13 CHL=OKN 
CKL=DKN 
RETURN I 

C 

END 
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SIEROUTINE HERAT (DK,RS) 
IiVPLICIT REALva(D) 



4 : 
* 
« 
* 
* 



THIS SUBROUTINE ITERATES THE FINITE DIFFERENCE 
EQUATION FROM THE BCTTON TO THE SURFACE. THE 
COMPUTED VARIABLES ARE: 



DUZ U) : VALUE OF THE MCCE EIGENFUNCTION 
MS : NUMBER OF MODE CROSSINGS 
* DUMAX = ABSOLUTE MAXIMUM VALUE OF EIGENFUNCTIO 



c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C<<<<<<<«<<<<<, <<<<<<<<<<<<<<<<<<<<<<<< COMMON ELCCKS >>»» 

C 

C 

CCMMCN /DBLl/ DVZ ( 501 )» DH, DH2 , DCBSG ,DRORB , CKAP , DUMAX 
CCMMON /DBLE2/ DW2 , DC B , COM I N , OW , C C I , DCl S Q 
COMMON /DBLE3/ DUZ ( 5 0 I ) , OZZ ( 50 I ) 

CCMMCN /PARAM/ N,NF1,IEX 

CCMMON /SINGl/ ZM , C B , CM 1 N , C MAX , P LT M AX , UM , F C 
CCMMON /DKMAP/ DKM I N , CKM AX , CKL , CKL 
CCMMON /DHHH/ DH2S IX , DH2D3 , CHHL F 
C 
C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DATA >>>>>>»>>>>>» 
C 

DATA DUPPER/0 .1060/ 

DATA NPRINT/6/ 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS »>» 

C 

C 

C DEFINE SCME CONSTANTS 

C 



0START=1.CC-10 
DK^— DK^DK 

CE=(DKAP-OK)*(DKAP+DK) 

C 

C r— SET UP RUN PARAMETERS 

C 

I DLMAX=OSTART 
DMARK= 1 .000 
MS = 0 
C 

C r INITIAL VALUES 

C 

DUJ=DSTART 

CLFJ=ORCRe*DSCRT ( ( DK-CKMI N ) ( DK + OKM IN ) )*DSTART 
DFJ=DE-DVZ(1) 

DLZm=CUJ 



C 

C-- — DEPTH LOOP BEGINS 

C 

c 

DC 8 J= 2 ,NP 1 
C 

C 'ITERATE THE FUNCTION 

C 

DFJP 1 =CE-DVZ( J ) 

DENOM= 1 . 0 DO+DH 2 SIX*CFJP 1 

OL JP 1 = ( ( 1 . 0 D 0 -DH 2 C 3 *CFJ )^-^CUJ + OH-CUPJ )/OENCM 
DA-( 1 . 0 DC-DH 2 C 3 ^DFJP 1 J-CUPJ 
CE=( LFJ-fDFJPl-DH 2 S IX=^ C F J-^OF JP U =»OUJ 
DLFJP 1 = (CA-D 6 ^DHHLF)/ 0 £NQM 
C 

C '-CHECK FOR MAX DU ( J ) 

C 
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CAUJPl = C/iBS(OUJPl) 

IF (OAUJPl-OUMAX) A, 4, 2 
2 IF (DAUJPI-LUPPER ) 3,3,5 
5 CLPAX=DAUJPl 
GC TO 7 

CHECK ZERC CROSSINGS 

A IF (DUJFl=i'DMARK) 5,6,7 

5 PS=KS+L 

CPARK=DSIGN(OMARK,CUJPl) 

GC TO 7 

6 IF (DUJPl) 7,5,7 

— ^ SET FOR NEXT ITERATION 

7 CFJ=DFJP1 
CLJ=DUJF1 
DLPJ=DUPJP1 
DL2( J)=DUJ 

8 CONTINUE 



DEPTH LOCF ENDS 



RETURN 



OVERFLOW PROTECT I CN 

REDUCE SCALE CF U(Z) 

9 IF (DSTART .LT.L.OD-60) GO TC 10 
0START=DSTART*1 .00-10 
GC TO 1 



10 WRITE (NPRINT,11I CK 
NS = C 
RETURN 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 
C 

11 FORMAT Cl SUBROUTINE ITERAT OVERFLOW PRCTECTICN*, 

^ • EXCEEDED. ', //, 

1 • MS SET TO ZERO AND RETURNED*, 

* • TO CALLING PROGRAM. ', //, 

2 • DK = • , G20. 12 ) 

END 
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SLEROUTINE FIXUP (DKN) 
IRFLICIT REAL«8(D) 



« THIS SUBROUTINE FORCES A CECAYING EXPONENTIAL 

^ ON THE EIGENFUNCTION ABOVE TFE UPPER TURNING 
^ POINT. 



<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CCMiMON BLCCKS >>>>>> 
CCRRON /PARAM/ N,NF1,IEX 

CCRMON /D3L1/ DV2 ( 50 1 ) , DH , DH2 t DC BSC ,OROR3 , C KA P , OUMAX 
CCRRON /CBLE3/ DU 2 ( 50 1 ) , DZZ ( 50 1 ) 

-CCMKON /DHHH/ DH2S IX, 0H203 ,DHHLF 

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FROGRAR BEGINS >»» 



DE=(DKAP-OKN)>!<(OKAP+DKN) ' 

NF2=NPl+l 
CAPG=DE-CVZ(NPL ) 

IF (DARG.GT.O.ODO) RETURN 
IFASS='L 

IF (DA3SIDUZ (NPl n .LT .DUMAX ) IPASS = i 
FIND THE UPPER TURNING POINT 



DC 1 1=2, NPl 
INDEX=NP2-I 
CKZ2=DE-CVZ( INDEX ) 

IF (DKZ2) 1,2,2 

1 CONTINUE 

RETURN 

^ DECAY TFE SCLUTION 

2 JTPU=INCEX 
JTPUP1=JTPU+1 
IST0P=NP1-JTPU 
CUJ=C.OCO 
CU'FJ=l.CD-5 • 

CFJ=DE-DV2(NP1) 

DC 3 J=1,IST0P 
INDEX=NP2-J 
CUZ( INCEX)=DUJ 
DFJP1=DE-DVZ( INDEX-1) 

DENOM= I .0D0+DH2S IX'^'DFJP I 

DUJP1 = ( (I .0D0-DH2C3=^CFJ )*CUJ + DH=s'CUPJ )/DENCR 
CA={ 1.0CC-0H2D5«DFJP1 )*DUPJ 
Ce=(CFJ+DFJF 1-DH2S I X^ CF J^DF JP 1 ) *DU J 
DUFJP1=(DA-C&*DHHLF)/CEN0M 
CFJ=DFJP1 
CLJ=0UJF1 
2 DUPJ=DUPJP1 

CCCEF = DUZ( JTPU/DUJPl 

DC A I=JTPUP1,N 
A OUZ( IJ = CUZ(I )^OCOEF 

IF (IPASS.GT.O) RETURN 
OUPAX=0.000 

CO 6 I=1,JTPU 
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C;iCU2I = CABS(0UZ( I ) ) 

IF (DACUZI-DUMAX) 6,6,5 

5 CLNAX=OADUZI 

6 CCMINUt 



RETURN 
C 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 



SteROUTINE CNORMl ( CIZ , CUMAX ,NP 1) 

IPFLICIT REAL’!^8(D) 

* 

♦ THIS SUBRCUTINE NGRNALI2ES ThE DOUBLE 

♦ PRECISION VECTOR CbZ (NPl), TC THE VALLE 
=«= DUMAX. 

♦ 

DINENSICN DUZ(NPl) 



C1EST=DLRAX*0.1D-60 



DC 1 1=1, NPl 
DLZI = DUZ1 I ) 

IF (DAdSIDUZI ) .LT.CTEST) DUZI=O.ODO 
1 DLZ(I)=DUZI/CUMAX 



RETURN 

END 
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SLEROUTINE BOTTOM (DK,MJ 
I^FLICIT REAL«S(D) 

c =* 

C => This SUBRCL'TINE COMPUTES A BOTTOM SCLND 

C * PROFILE FOR PLOTTING PURPOSES. 

C =^‘ 

Q A is 4: *>>;****&#;!; 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLOCKS >»>» 
C 

CCMMCN /DBLl/ DV Z ( 50H , CH , CF2 , DC 6 S C , DROR8 t C K A P , DU M AX 
COMMON /DBLE2/ OW2 , DCB , OCM I N ,Dl’< ,CCi ,OCl SC 
COMMON /DELE3/ DU Z 1 50 L ) , DZZ ( 50 1 ) 

COMMON /SINGl/ ZM , C B , C M I N , C MAX , F LT M AX , U M , FC 
CCMKON /BOTTO/ DEC AY , ZB ( 20 ) , UB ( 2C ) 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM EEGINS >»» 
C 

C COMPUTE CEPTh VALUES 

C 

C 

FE= (ZM-PLT MAX 1*0 .5E-1 
DF.B=OBLE(FB) 

C 

C COMPUTE SOME CONSTANTS 

C 

CGAM=DSQRT( ( DK-D’»^/DCE J »!■ ( DK+DW/DC B )) 

DARG=DGAM--DHB 
IF CDEXFU+OARG) 3,3,L 

1 CEX=CEXP(DARG) 

CV = DUZ( 1 J'-^-'CRCRB 

C 

C : COMPUTE THE PROFILE 

C 

C 

DC 2 1=1, 2C 

IF (DV.LT.l.ODO-30) GC TO 4 
DV=CV^DEX 

2 UE(I)=SNGL(DV ) 

C 

RETURN 

C 

C-^ . NULL PROFILE 

C 

3 1=1 
C 

A CC 5 J=I,2C 
5 UE(JI=0.0 
C 

RETURN 

END 
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SLERCUTINE INTEGR (CK,RCOE) 

IRFLICIT REAL*8(0) 

^ >5: 5?: ^ is 3{C 3?: 351 ^ a?f ;i« A ^ xi 5ft ^ V * 3^- V* 3j» 4; ^ r; 

c * ' ' 

C * THIS SUBROUTINE COMPUTES THE NC RM AL I Z AT I C N 

C * INTEGRAL, AND THE SOUND VELOCITY NORMAL- 

C * IZATICN INTEGRAL. THEN THE GRCUP AND PHASE 

C VELOCITIES ARE COMPUTED. 

C * 

c 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLOCKS >>>>>> 
C 

c 

COMMON /PARAM/ N,NP1,IEX 

COMMON /SINGl/ ZM , CB , CM I N , C M AX , PLT MAX ,UM , F 0 
CCMiMON /DELI/ OV^Z ( 50 L ), LH , DH 2 , OCdSO ,DRORE , C KAP , DUMAX 
COMMON /DBLE2/ DW2 , CC B , COM I N , DW , CC 1 , DCl SC 
COMMON /DBLE3/ DUZ ( 501 ) ,OZZ (501 ) 

COMMON /DELE4/ D S , C SSCC , DP V EL , DG V EL 
COMMON /RHC/ RO,RB 
COMMON /OSCUNO/ DCC(501) 

COMMON /BCTTO/ DECAY, ZB ( 20 ) , U6 ( 2C ) 

COMMON /CKMAP/ OKM I N , CKM AX , CKL , OKU 
ECUIVALENCE (OUZl,CUZ(in 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DATA >>>>>>>>>>>>>> 
C 

DATA NREAC,NPRINT/5,6/ 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS >»» 
C 

c 

C *** SET UP THE END POINTS FOR THE TRAPA2CIDAL RULE 

C 

DLZ1SG=CUZ(1 ) 

IF (DABS(OUZISG) .LT. l.CC-30) OUZISG=O.ODO 

CUZ1SQ=CU21SC^DUZISC 

CS=C.5CCv(DUZLSD) 

DCISQ=CCC( 1) *OCC( I ) 

CUZ1SG=CUZ1SQ/(DC1SG) 

CSSCC=0.500*(DUZ1SG) 

c 

Q INTEGRATING LOOP BEGINS 

C 

c 

CC 1 1=2, N 
CL2I2=CLZ( I) 

IF (DAeS(DUZI2) .LT.l. CD-30) GO TC I 

CLZI2=CUZI2*CUZI2 

CS=DS+0LZI2 

CCCI=DCC(I) 




C 

C 

C INTEGRATING LOOP ENDS 

C 

C 

C *=«■* ADD BOTTOM EXPONENTIAL DECAY *** 

C 

CA2=( OK-CKMIN ( DK + CKM I N ) 

CE = CBLE (RB)^OUZl vC.5DC/0SQRT(DA2 ) 

CS = CELE (RC )*DH=«=OS + Ce 
CSSCC=DSSCC*DBLE( RC)*CH+D6/CC6SG 
DECAY=CW=*'C6/( DCB^^OK^OS ) 

C 

C 

C COMPUTE THE PHASE AND GROUP VELOCITIES 

C 
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CFVEL=Cl%/DK 

DGVEL = CS/ (CP VEIL'D S SCO 
C 

wPITE (NPRINT,2) PGCE.OPVEL 
UPITE (NPRINT,3) PGDE.DGVEL 
C 
C 

RETURN 

C 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 



2 FORMAT 


( • OPHASE 


VELCCITY 


FOR 


MOCE* , 


14, 


' IS ', 


1 


G15.7, • 


M/SEC. *J 










2 FORMAT 


{ 'OC-RGUP 


VELCCITY 


FOR 


MOCE • » 


14, 


' IS ', 



1 G15.7, • P/SEC.') 

ENC 
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SieROUTINE MODPLT (M,JUMP) 
IPFLICIT REAL^8(D) 



« 

* THIS SUBROUTINE PLOTS THE MODE SHAPE 

=«« CN THE PRINTER USING THE N.P.S. ROUTINE 

* UTPLOT . 

* 

=« IF JUMP=1, ALL POINTS IN THE MODE 

* FUNCTION VvILL EE PLCTTED. IF JUHF 

’i' IS NOT EQUAL TO 1, EVERY NPl/lOO TH 

^ POINT WILL BE PLOTTED. 

ft :jc ic ft :(( >c ^ ;;; >i: ^ 



COMMON BLOCKS »>>»> 

CCMMCN /D6L1/ DV Z ( 50 1 J , CH , DH2 , DC E SC , CRORB , C K AP2 , DUP AX 
CCPMON /DBLE2/ DW 2 , CC B t CCM I N, DW t DC 1 ,0C L SQ 
CCPPCN /DBLEB/ DUZ ( 5C I ) , DZZ (50 L ) 

CCMHON /PARAM/ N,NP1,I£X 
CCPMON /HEAD/ HED(20) 

CCPMON /SINGl/ ZM,Ce,CMIN,CPAX,PLTPAX,UM,FG 
COMMON /DMAP/ DKl(lOO) 

COMMON /FL0T6/ RANGE(^) 

COMMON /INPUT/ NUMV,NMCD 

COMMON /BCTTO/ D EC A Y , ZB ( 20 ) , UB ( 2 C ) 

PROGRAM BEGINS >>>>>> 



CK=DK1( M ) 

CHECK IF FULL PLOT OF ALL POINTS DESIFEC *=** 

IF (JUMP-1) 1,2,1 

1 JLMP=NP1/1G0 

2 JLMP2 = JUMP’!>2 
NLMB=NP1/JUMP 

PRINT THE HEACINGS ** 

WRITE (6,6) HED,FG 
WRITE (6,7) M,DKl(M) 

WRITE (6,8) 

ftftft COMPUTE THE RANGE FACTORS FOP PLOT SIZE 

PANGE(1)=1.1 
RANGE(2)=-l.l 
RANGEO )=2M 

WILL BOTTOM PROFILE BE INCLUDED? 

IF (ZM-PLTMAX) 4,3,3 

3 M0CCUR=O 
RANGE(4 ) = C.O 
GO TO 5 

draw bottom PROFILE 

4 MCCCUR=1 
RANG£(4)=ZM-PLTMAX 

CALL UTPLOT ( UB , ZE , 20 , RANGE , 1 , MODCUR ) 

M0DCUR=3 



DRAW OCEAN U(Z) PROFILE 
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c 

5 C/SLL UTPLQT ( CUZ , CZ2 , NUMB , RANGE , JUMP2 , MOCCLR ) 

C 

C 

RETURN 

C 

6 FCRMAT CP, 20A4, • FREQUENCY : F6.1, ' HZ') 

7 FCRRAT CO', T25, 'MODE', 14, ' : K =', F16.12) 

8 FCRMAT C.O', // , T23, 'NQRMALIZEC E I GEN FUNC T I CN 

I • VS DEPTH') 

END 
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SL6R0UTINE OUTPUT ( CU Z , CK , MODE , N FTS ) 

IMPLICIT REAL^Sin) 

CIMENSICN CZZ(NPTS), DUZ(NPTS), CCC(NPTS) 

CCMMON /RHC/ RQ,RE 

CCMMON /SINGl/ ZM , C8 , CM I N , CMAX , P LT MAX , U M , F C 
CCMMGN /BCTTO/ DECAY, ZB ( 20 ) ,UB ( 20 ) 

CCMMQN /DBLE4/ D S , D SSCC , OP V E L , DC V E L 
CCMMON /HEAD/ HED(20) 

CCMMQN /WCRK/ EXTRA(20) 

DATA NPRINT,NPUNCH/6,7/ 

C 

C 

1 kRITE (NPUNCHtlO) MODE , FQ , DK ,DPVE L , OGVEL , CS , DEC AY 
C 

DC 2 1=1,20 

2 EXTRA! I )=UB(21-I ) 

C 

WRITE (KPUNCH,m EXTRA 
WRITE (NPUNCH,U) CUZ 
WRITE (NPRINT,13) MQDE 
RETURN 
C 
C 

ENTRY PROFILI CZZ,CCC,NPTS ) 

WRITE (NPUNCH,?) HED 
NFTS20=NPTS+20 

WRITE (NPUNCH,8) N FTS20 , RO , R E, CM I N , CB , ZM , PLTMAX 
C 

DC 3 1=1,20 

3 EXTRAd ) = ZM-ZB(21-I ) 

C 

WRITE (NPUNCH,12I EXTRA 
C2M=06LE(ZM) 

C 

DC 4 I=l,NPTS 

4 C22( I)=D2M'*DZZ( I ) 

C 

WRITE (NPUNICH,12) DZZ 
C 

DC 5 1=1,20 

5 EXTRA(I)=CR 
C 

WRITE (NPUNICH,12) EXTRA 
WRITE (NPUNCH,12) CCC 
WRITE (NPRINT,9) 

C 

DC 6 I=1,NFTS 

6 DZZ( I) = DZM'OZZ( I ) 

C 

RETURN 

C 

7 FORMAT (20A4) 

e FCRMAT (I10,2F10.5,4F10.2) 

9 FCRMAT (‘OSCUND VELOCITY DETAILED PROFILE PUNCHED* ) 
END 
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SLbRCUTIKE CZFLOT 
IMPLICIT REAL=4=8(D) 

C 

C « THIS rSUBRGUTINE PLCTS THE SCLND VELCCITV 

C * PROFILE CN THE PRINTER PLOT USING 

C * The N.P.S. ROUTINE UTPLCT. 

C * 

Q + j{! jJ; * :ic :(e :j« ;{0{< * * A *★* 5)t ^ ^ 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLCCKS >>>» 
C 

CCMMON /DBLl/ DVZ ( 50 1 ) , OH , DH2 , DCBS C ,DRORB , C K A F2 , DUN AX 
CCPMON /0eLE2/ DK2 , CC B , DCM I N , D W , DC 1 ,OC 1 S G 
CCNMON /DBLE3/ DUZ ( 5 C 1 ) »DZZ ( 50 I ) 

COMMON /PARAM/ N,NP1,IEX 
COMMON /HEAD/ HEC(20J 
COMMON /UCRK/ WRKU20J 

COMMON /SINGL/ ZM , 0 B , CM IN , CMAX , P LTMAX, UM , FG 
COMMON /PLCT6/ RANGEI4) 

COMMON /DSCUNO/ 0CC(50U 

COMMON /EGTTO/ DECAY, ZB ( 20 ) ,UB ( 20 ) 

C 

C«<<<<<<<<<<<<<««<<<<<«<«<««<<< PROGRAM BEGINS >»»> 
C 

C *** URITE THE HEADINGS 

C 

k\RITE ( 6 , 5 ) HED 
^RITE (6,6) 

WRITE (6,7) 

C 

IF (PLTMAX.lt. ZM) PLTMAX=ZM 
FE=(ZM-FLTMAXJ“0.5C-1 

c 

c — DEPTH INCREMENT LCCP BEGINS 

C 

DC 1 1=1,20 
WPK 1( I )=CB 
ZE(i)=HE«fLOAT( I-l ) 

1 CONTINUE 
C 

Q DEPTH LGCF ENDS 

C 

C =!=## COMPUTE THE RANGE FACTORS 

' C 

RANGE! 1 ) = AMAXl(CB,CMAX) 

RANGE(2 ) = CMIN 
RANGEI3 )=2M- 
-C 

C *** WILL BOTTOM PROFILE BE DRAWN? 

C 

IF (ZM-PLTMAX) 3,2,2 
C 
C 

2 MCDCUR=0 
RANGE(4)=0.0 
GC TO 4 

C 

3 MCCCUR=1 

RANGE(4 J = ZM-PLTMAX 
C 

C DRAW THE BOTTOM C(Z) PROFILE 

C 

CALL UTPLOT ( WRK 1 , ZB , 20 , RAN GE , 1 , MCDCUR ) 

MCCCUR=3 

C 

C DRAW THE OCEAN C(Z) PROFILE *=** 

C 

C 

4 CALL UTPLCT ( CCC , CZZ , NP 1 , RANGE , 2 , MCDCUR ) 

C 

RETURN 
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5 FCPNAT (M‘, 20A4) 

6 FCRMAT CC'f // , T28, 'PLOT CF C(Z) \/S CEFTHM 

7 FCRMAT COS //, Til, 'BOTTOM AT ZERO. DEPTH IN 
1 'METERS AECVE BOTTOM. C(Z) IN METERS/SEC) 

ENC 
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I^PLICIT REALMS (A-H, V-Z ) , L0GICAL*4 ($) 

C 

c 

c 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<COMMOK BLOCKS >>>>>>>>>>>>>>> 
C 

VEL0(1021), DEPTH (1021) 

ZZZ( 1021) , DVZ( lOCl) 

CC(50) ,DEP(50 ) 

DRn,DRB,DRORB 
CMAX,CMIK,CB,CBSQ 
H, CHfHPLOT ,OhhLF, DH2SIX,0H2C3, 

DH2,H0OTT ,DHB 

DKtCKMIK ,DKMAX ,DKOLC,CKU, DKL, 
DTEST,DSMAX,0£ 

ZMAX,CSS,DECAY,OPVEL,DGVEL tGAM 
WOCe,WCCNPl,WOCMIN,CVZB 
FCV (20) ,FG ,W, ^^2 
DKM,CINT,DIFF,CS, DM 
DKl ( 100) 

H£D(2C) 

K‘JMV,NMOO,NFREC,KMCCF,N,NP1, 
NP21,KFTS,MCD£,IT ,KF2 
LOMODE, IHMQDE ,KUELL , I WELL , KV. E LL » 

JHI , JLCW 

JUPPER , JBCTT , ISHAF , JTPL, JTPU 
NREAD,NPRINT,NPUNCF. 

$GRAPF, a.BCTPR,TCARCS,$CEL 
DPI ,DPIC2,DPIC4,DPIRT2,D2PI 
C 

C ' 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS >>>>>>>>>>>>>» 

C 

C 

WRITE (KPRINT,6) 

CALL INFUT2 
CALL LINEAR 

IF ($GRAPH) CALL CZPLCT 

IF ($CEL.OR,$CAROS) CALL PRCFIL ( C EPTH , VE LC , NPTS ) 



CCMMCN 


/ 


D1 


/ 


COMMON 


/ 


02 


/ 


COMMON 


/ 


D3 


/ 


COMMON 


/ 


D4 


/ 


COMMON 


/ 


C5 


/ 


COMMON 


/ 


06 


/ 


COMMON 


/ 


D7 


/ 


COMMON 


/ 


OS 


/ 


COMMON 


/ 


010/ 


COMMON 


/ 


DW 


/ 


COMMON 


/AWKB/ 


CCMMON/DKMAP/ 


COMMON 


/HEAD/ 


COMMON 


/ 


11 


/ 


COMMON 


/ 


•12 


/ 


COMMON 


/ 


13 


/ 


COMMON 


/ 


10 


/ 


COMMON 


/ 


LI 


/ 


COMMON 


/ 


PI 


/ 



C 

c- 

c 



DC 5 IFREC=1,NFREC 



•FREQUENCY LCCP BEGINS 



C 

C- 

C 



FC = FQV( IFREQ) 

W=FC*D2FI 

W2=W«W 

wcce=w/ce 

WCCMIN=W/CMIN 



CC 1 J=l,NPl 
UCC=W/VELC ( J ) 

1 DVZ( J)=(WCCMIN-WCC)*( VCC+WOCMIN) 
WCCNPl = VsCC 

cvze=( WCCM IN-WOCB I^^'IWCCMIN+WOCB) 
CALL MAXMI3 165) 



•COMPUTE VELOCITY FUNCTION 
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no oooo no o ooo 



KWELL=l 

LCNCDE=C 

^QQE LOOP BEGINS 

CC A M0DE=1,NM00F 
C/5LL SE/JRCh 
C/iLL ShAPE4 
CALL INTEGR 
IF (JGRAPhI GO TO 2 

IF (HCCE.EC.U rtRITE (NPRINT,?) FED 

IvRITE (KPRINTtS) R'ODE , CK , DP VEL , OG VEL , DS S , C ECAY 

GC TO 3 

2 CALL MODPLT 

3 IF ($CEL.OR.$CARDS) CALL OUTPUT (ZZZ,NPTSJ 
A CONTINUE 

5 CONTINUE 



FREQUENCY AND MODE LCCFS END 

STCP 



C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 
C 

6 FORMAT Cl N0RM04 : WKB METHOD NORMAL MODE', 

1 • EIGENVALUE SOLVER. • ) 

7 FORMAT ( • I* , ICA8 ,//, • MODE • , T2 5 , • FOR I Z RAVE N0.',T40, 
I'FHASE SP EEC* ,T60, 'GROUP SP E ED • , T E4 , • DSS • , T 102 , 

8 FORMAT ('O' ,I4,F20.12,' 1/M ET ERS • ,T40 , F7 . 1 , • M/SEC, 
1T60,F7.1,‘ M/SEC ,T76,2G20. 12) 

ENC 
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ELCCK CATA 

iNPLICIT BEAL«3 (A-H, V-Z ) , L0GICAL*4 {$) 
CCMMUN / D5 / CMAX,CMIN,Ca ,CI3SQ 
CCPMQN /AWKB/ DK.\ , D INT , C I F r , DS , DN 
CC^MON / PI / DPI ,DP1C2 ,DPIC4tDPIPT2,02PI 

CCMMON / II / NUMV,NMOD,NFREC,N!^CDF,N,NPlf 

3 NP21 ,NPTS ,MOCE, IT ,NP2 

CCMMON / IG / NREAD,NPRINT,NPUNCF 

CCMHON / LI / $GRAPH, $BCTPR,$CARCS,SCEL 

DATA NREAC, NPRINT, NPUNCH /5, 6, 7/ 

DATA $oRAPH,$CAROS,$CEL,$BCTPR /3*. FALSE., 
DATA CHAXtCMIN /O.CCC, 9.9D9/ 

DATA MODE /O/, OP /C.OCO/ 

DATA 0PI02 /I .5707962267948966/ 

DATA D2FI /6 .283 1852C7 179586/ 

DATA DPI04 /0 .785398163397^483/ 

DATA DPI /2. 141592653589793/ 

END 



.TRUE./ 
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I 



It 



no o no on o ooo ooo oooonoo 



SLeRQUTINE KAXMI3 ( ’> ) 

UPLICIT REAL=<=8 (A-h, V-Z ) , LOGICAL*^ ($) 



LAST CHANGED 3 AUG 73 



COMMON BLCCKS >>>>>> 



C 

C 



COMMON / C5 / CMAX,CHIN,C6,C3SQ 
COMMON / D7 / OK,CKMlN,DKMAX,OKCLC,CKUf DKLt 
; DTEST,DSM;AX,DE 

COMMON / CV, / FCV(20J ,FC,U,U2 
COMMON /AWKB/ DKN , D I NT , C I FF , OS , C M 
COMMON /HEAD/ HEC(20) 

COMMON / II / NUMV,NMOC,NFRECtNMCCF,N,NPl, 

5 NP21 tNPTS,MOD£, I T ,NF2 

COMMON / 10 / NREAC,NPRINT,NPUNCF 
COMMON / LI / $GRAPr, $BCTPR tiCARCS ,$CEL 
COMMON / PI / DPI , DPI02,DPIC4,DFIRT2,D2PI 



PROGRAM BEGINS >>>>> 



CM=O.OCO 

OKMIN=W/CB 

CKMAX=W/CMIN 

CKN=OKMIN 

CALL WKE 

CSMAX=DS/OPI 

MCCMAX=ICIi\T(DSMAX) 

CKCLD=DKMAX 

WRITE (NPRINT.3) HED.FQ 

CKN=DKMAX 

MCDMIN=C 

WRITE (NPRINT,4) 

DATA MAXI/'MAXI '/, MI M/ 'MINI*/ 
WRITE (NPRINT.5) M AXI , DKMAX , MCOM I N 
WRITE (NPRINT,5) M I N I ♦ DKM I N , MOD M AX 



NMCDF = i\MCC 

IF (MODMAX-NMOO) Ii2,2 



•CHECK NUMBER MODES RECUESTED 



1 CONTINUE 



■IF WE GET FERE WE MOST RESET 
THE NUMBER CF MODES RECUESTED 



NMCDF=MCDMAX 

WRITE (NPRINT,6) 
WRITE (NPRINT,7) 
WRITE (NPRINT,8) 
WRITE (NPRINT,7) 
WRITE (NPRINT,6) 
WRITE (NPRINT.9) 



NMCDF 



2 IF (MOCMAX.NE.O) RETURN 



•IF wE GET FERE WE ARE BELOW 
CUT - OFF FREQUENCY 



WRITE (NPPINT,6) 
WRITE (NPRINT,7) 
WRITE (NPRINTtO) 
WRITE (NPRINT,7) 
WRITE (NPR1NT,6) 
WRITE (NPRINT.IO) 

RETURN I 



135 






c 

c 



FCRMAT Cl', lOAe, AX, 'FRECUENCY : F6.1J 

FCFMAK 'C ,// ,T2£>, '6CUNCS ON HCRIZCNTAL ^.AVE NUMBER') 
FCRMAT CC, T22, AA, 'MUM K :', F 12 . 8 , 3 X , ' MCOE S :',IA 
FCRMAT (T33, 13( '* ' ) ) 

FCRMAT (T33, IF^, 11X,IH=») 

FCRMAT (T32,'* CALTIGN * • ) 

FCRMAT ('O', T26, 'MODES REQUESTEC RESET TG lA) 

1C FCRMAT ('0',/,T25, 'FRECUENCY RECUESTEO IS EELOW , 

1 ' CUT-OFF ',// ,T20, ' JOB TERMINATED OR NEVx FREQ', 

2 'LENCY SELECTED' ) 

END 
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SL6RCUTINE IMEGR 

I^PLICIT REALMS (A-H, V-Z), L0GICAL*4 ($) 

LAST CHANGED 9 AUG 73 



CCf'l^QN / 
CCNiMON / 
CCP'^O.M / 
CCPMGN / 
CCR'MOK / 

1 

CCPPON / 

2 

CCMMON / 
CCRPON / 
CCP/''ON / 
CCMMON / 
3 



D1 / VEL0(1021), DEPTH (1021J 
C2 / ZZZ( 1C21 ) , DVZ( lUOl ) 

D4 / 0R0,CR&,0R0R3 
C5 / CMAX ,C-y IN,CB, cesQ 
06 / H, CH thPLCT ,DHhLF , DH2SIX, DH2C3 , 
DH2,HE0TT,0HB 

C7 / OK ,CKKIK,CKMAX ,CK0LC,CKU, DKL , 
DTEST,OSMAX,DE 

D9 / ZMAXfCSSf DECAY, DPVELtDGVEL, GAN 
010/ WOCe ,l^CC^Pl ,WCCNIN ,CVZb 
CW / FCV (20) ,FC,W, W2 
II / NUMV,i\MOD,NFREC,NNCDF,N,NPl , 
NP2l,i\PTS,M0DE, IT ,NP2 



DSS=C.0C0 

C.SSCC=0.000 



<<<<<<<<<<<<<<<<<<<<<<<<<<<<< INTEGRATING LOOP BEGINS 

CC I 1=2, N 
ZZZI2=ZZZ( I ) 

IF (0AES1Z2ZI2) .LT. 1.00-30) GO TO 1 
ZZZI2 = ZZZ12="-ZZZI2 

CSS=DSS+ZZZI2 

VELCI = VELC( I ) 

2ZZI2 = ZZZI2/(VELQI*V£LCI ) 

CSSCC=DSSCC+ZZZI2 

1 CONTINUE 

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< INTEGRATING LCCP ENOS 



IF (DABS (ZZZ(NP1 )) .LT .1 .00-30) GO TO 3 
A2=DSQRT(CVZ3-DE) 

e = ZZZ(NPl ) vZZZ(NPL)’*=DR0R6^DR0/ ( 2 .0D0-A2 ) 

2 CSS=CRO=^''Oh^"-CSS+B 
CSSCC=DSSCC='^DRO = OH + B/CaSQ 
C£CAY=W^fc/ (C8*DK*DSS ) 

C 

CPVEL=W/OK 

CGVEL = OSS/ (OPV£L=«CSSCC) 

RETURN 

C 

3 E=C.000 
GO TO 2 

C 

END 
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SIBROUTINE WKB 

I^PLIC^T PEAL-8 ( A-H, V-Z) fLCGICAL'-i-H ($) 

C 

C 

LAST CHANGED : 7 SEPT 73 ASY.MPTCTIC B. C. 



CCMMQN ELCCKS >>»» 



CCPPON / D2 / 
CCPMON / DA / 
CCPMUN / D6 / 

1 

CCMMQN / 07 / 

2 

CCMMQN / 010/ 
CCMMQN /Al^Kd/ 
CCMMQN / II / 

3 

CCMMQN / 12 / 

4 

CCMMQN / 13 / 
CCMMQN / LI / 
COMMON / PI / 
COMMON / OTH/ 
CATA DEXPU /I 



ZZZ(1021), DVZ(LOOl) 
ORGfCRBi-ORORB 

H, CH,HPLCT,OHELF, OH2SIX,OH2C3 , 
0H2 ,HB0TT ,0H8 

OK ,DKMIN ,0KMAX rOKQLC ,0KU ,0KL i 

OTESTf OSMAX.OE 

WOCB ,V\CCMPL fWOCMI N ,CVZB 

0KN,DINT,0IFF ,C5,0M 

NUMV , NMQC , NFRtC,NMOCF,N,NPl, 

NP21,f\PTS,MGDE,lT,NF2 

LOMIOCE, IHM00E,KWELL tlivELL i MW ELL, 

JHI , JLOW 

JUPPER, jaCTT , ISHAP , JTPL, JTFL 
$GRAPF, $aQTPR,$CARCS,$CEL 
DPI ,0FI02,OPI04,DFIPT2,02PI 
DTHETL 

.702/, DWELL /6.0/ 



PROGRAM BEGINS »>>> 



CINT=O.COO 

JTFL=l 

CE= (DKMAX-OKM* (DKMAX + CKN) 

CKB=DSQRT( DVZB-OE) 

$TOP=. FALSE. 

IWELL=l 

ISIG=-l 

CAPG=O.OCO 

CARGP=C.CDC 

CEPTH LOOP BEGINS 



CO 21 J=1,NP1 
CKZ2=D£-DVZ( J ) 



1 

2 
3 



IF 


(DKZ2) 


1,2,3 


IF 


(ISIG) 


4,5,6 


IF 


( ISIG) 


7,21,9 


IF 


(ISIG) 


10,19,20 



CHECK FOR CONDITIONS 
AND TYPE OF SOLUTION 



IMAGINARY REGION 



CARGP=DSQRT(-DKZ2J 
CAPG=CARG+CARGP 
GO TO 21 

5 ISIG=-1 
GO TO 4 

6 CARGP=DSGRT(-DKZ2) 

DX=CKZ/ (CKZ+OARGP) 

C INT=DIM.T.+ OKZ^(DX-i. 000)- 0. SCOTCH 

IS IG=-1 

JTFL^J-1 

CARG={ 2.00 C-DX)=^DARGP* C. 500 
GC TO 21 



ZERO VERTICAL t^AVE NUMBER 
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7 IE/1CK = -1 

8 ISIG=0 
GC TO 

S IE/iCK=l 
GC TO 8 

* peal region 

10 ISIG=l 

IF ($TCF) GC TO 11 
JTPL=J 

11 CKZ=DScRT(DKZ2) 

CX=CARGF/ (LKZ+DARGF ) 

C>5RG=DARG-*DARGP«{ DX-l .CC0)*C.5DC 
IF ($T0P) GO TO 13 

12 CThETU=CATAMDTANF(CARG=i'CH) ) 

C INT = DTFETU 

$TCP= .TRUE. 

GC TO IE 

13 CT l = DMOD( 0INT,D2PI ) 

IF (DARG*Dh-DWELL ) 15tlA,14 

IF (IWELL.EG.KWELL) GO TO 28 

C/!PG=C.CC0 

IV>ELL=IUELL + 1 

JTPU=J 

CTFETU=CPIQ4 
C INT=DTFETU 
GC TO IE 
15 A = CSI.NI(DT11 
e=CC0S( CTl ) 

CEL = DEXP(DARG-^DH) 

C£PL = DEXP(-DARG--Ch) 

AC=CEL-(A+E )+CEML* (A-e ) 

BC = DEL=MA + 8)-DEML*( A-B) 

D INTP=DINT-OT 1 + DAT AN 2 ( AC , BO } 

U IF (DINTP.GE.DIMT-CTESTI GO TO 17 
CINTP=DINTP+DPi02 
GC TO 16 

17 DINT=DINTP 

18 CINT = 0INT+t2.C00^CX)^DKZ*C.5D0’i^DH 
CARG=0.CC0 

GC TO 21 

19 CKZ=DSCRT (CKZ2) 

I SIG= + 1 

C INT=DINT+OKZ^DH 
IF ($TGF) GO TO 13 
JTPL‘=J 
GC TO 12 

2C CKZ=DSCRT(DKZ2) 

CINT=DINT+OKZ’^DH 

21 CONTINUE 



DEPTH LCCF ENDS 

IF (DKZ2) 22,27,25 

DECAYING SCLUTION 

22 CKS=DARGP 

DA RG= (DARG-C .500^0 ARC-P )*DH*2.0D0 
IF (DARC— DEXPU) 23 ,24,2^* 

2 3 C2=( DKS-nR0R6“-DNB} *CEXP (-DARG)/ { CKS + DR0RE*CKB) 
CTHETL=CATAN( ( 1.0CC+D2)/( 1.0DO-D2) ) 

GC TO 27 
24 CTFETL=CPIQ4 
GC TO 27 
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C -BOTTOM 

c 

25 JTFL=NPl 

IF (DKB.6C.C.0D0) GO TO 26 
CThETL = CATAN(DKZ/ { CKB’i'CRORB ) ) 

GC TO 27 
C 

26 CThETL=CFIC2 



COMPUTE 



27 CS = CINT + DTI-ETL 
CINT=DINT-DTHETU 
C IFF = DS-DM-^DPI 
^^^ELL = IWELL 
RETURN 
C 

26 IkELL=IWELL+l 
•GC TO 2A 
C 

ENC 



REFLECTED 



CIFFERENCE 
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SLEROUTINE SEARCH 

IMPLICIT ftEAL ’8 (A-H, V-Z), L0GICAL*4 ($) 

C 

C — — — — 

C LAST CHANGE : 15 AUGUST 1973 

Q 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON BLOCKS >>»» 
C 



CCMMON 


/ 


07 


/ 


OK, DKMIN fDKHAX,DKCLC , OKU, OKL , 


2 








OTEST,DSMAX,DE 


CCMMON 


/ 


06 


/ 


OIFF I 


CCMMON/ CKMAP/ 


OK 1 ( 100 ) 


CCMMiON 


/A'aKB/ 


OKN,CINT,CIFF,CS, CM 


CCMMON 


/ 


I 1 


/ 


NUMV ,NMCD ,NFREG,NMCCF ,N,NPl , 


3 








NP21,NPTS,M0DE, IT ,NP2 


CCMMON 


/ 


•12 


/ 


LOMGCE , IHMGDE ,KWELL , IWELL ,NWELL, 


4 








JHI , JLQ'W 


CCMMON 


/ 


13 


/ 


JUPPER, JBGTT, ISHAPtJTPL, JTPU 


CCMMON 


/ 


in 


/ 


NREAOtNPRINT ,i\PUNCH 


CCMMON 


/ 


L 1 


/ 


SGRAPH, $eCTPR, SCARCS i $CEL 


CCMMON 


/ 


PI 


/ 


DPI , CFI02,DPIC4, DFIBT2, 02 PI 


CCMMON 


/ 


OTH/ 


DTHETU 


CATA ITMAX , 


/30/ 



C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM BEGINS >>>» 

C 

C 

C . INITIAL SETTINGS 

C 

FIX=1.4DC 

IF (MOCE.GT.l) GO TO 1 

KkELL=l 

$LCW=. FALSE. 

$FI = . FALSE . 

CKCLD=DKHAX 
JLPPER= 1 
JFI = 1 
JLCW=NP1 
IFMCDE=C 
LCMCDE=C 
CKCEN=CKOLD 
CKC£N2=CKCEN 
JBCTT=NP1 
C 

C SET UP FOR FIRST GUESS 

C 

1 CM=DFLOAT(MODE-LOMODE-IHMOOE) 

IT = 0 

• IF ^$LOl^) GO TO 11 

CKGUES= (CKGLD-OKM IN ) / ( DSMAX-DM+ 1 .000) 

2 CKL=CKCLC-CKGUES*FIX 

IF (OKL .LT .DKMIN) CKL=CKMIN 
CXU=CKQLD 
FL=-DPI 
2 CKN=OKL 
CALL WKB 

IF (JTPL.LT.JhI) GC TC 5 
IF (OIFF.LT.O.ODO) GO TO 4 
FL=CIFF 
JHI=JTPL 
JLCH=JTFL 
GC TO 16 
C 
C 

4 F IX = FIX + C .400 
GC TO 2 
C 
C 

C — UPPER CHANNEL SET-UP 

C 



141 



ooo o o noo ooo non ooo 



c 



5 IF {DS-CFLCAT{ IHMCCE + l )*DPI ) 6,32,32 
t KInELL = 2 
6C TO 3 



7 

e 






1C 



IF I=. TRUE. 

JIFFER=JTPL 
IFRCOE=IhPCDE+l 
CR = DF10AT( IhMOOE ) 
CKN=(CKU+DKLJ*O.5D0 
CALL WKB 

IF ( JTPL.LT. JUPPER ) GO TO 9 

CKU=0KN 

GC TO 7 

IF (DIFF.LT.C.ODC) GO TO 10 

CKN=( CKU+0KN)’!'0.5C0 

GC TO 8 

CKU=DKN 

FL=CIFF 

GC TO 16 



LOWER CHANNEL SET-LP 



11 CKL=DK 
CN=CFLCAT (LCMODE) 

CI<L = DKCEN2 

CKN = {DNU + DKL)=?=0.50C 

12 CALL WKB 

IF (JTPL.GT.JBOTT) GO TC 13 
DKN=( DKt+DKiM)'‘0.5DC 
GC TO 12 
12 FL=DIFF 
CKU=DKI\1 

lA IF (FL .GT.O.OCO) GC TO 16 

CKL=DKL-.21:0C-(DKRAX-DKMN)/DSMAX 
IF (OKL .LT .DKMIN) CKL=0KMIN 
CI<N=DKL 
CALL WKB 

IF (JTPU.LT.JLOW) GO TO 15 

FL=CIFF 

GC TO lA 

15 KWELL=i 
$LCW=. FALSE. 

LCyCDE=LGMCDE-l 
GC TO 1 

CCRPUTE NEXT DKN 

16 DKN=DKL*FL'‘(CKU-CKU/ (FL-FU ) 
^CALL WKB 

17 CALL WKB 
IT=IT+L 

IF ($LOW.OR.$HI ) GC TC 18 

IF (JTPU. LT. JUPPER ) C I FF = D I FF-D F LGAT ( I HMCC E ) *CP I 
IF (JTPL.GT.JBOTT) D I FF=D I FF-DF LO AT ( LOMOC E J =»CP I 

CHECK WHETHER ABOVE CR BELOW 

16 IF (DIFFJ 19,20,21 
19 CKL=DKM 
FL=DIFF 
GC TO 22 

2C WRITE (NPRINT,31) PCDE 
GC TO 25 

21 CKL=DKN 
f L=DIFF 



TEST THE RESULT 



1A2 



ooo non non non non non 



22 CINC = 0/i3S(CKGLD-DKN)/DKN 
DKCLD=DKN 

IF (DAfaS(CirF)-DTEST) 25,25,23 

23 IF (M0U(IT,10).EC.9) GQ TO 2^ 

IF ( IT.LE.ITHAX) GC TC 16 

IF (DINC.LT. I .00-15) GG TO 25 
2k CKN = ( CKL + 0KL)’^0.5C0 
GC TO 17 

— 71-E ii^cDE HAS BEEN FCUNC 

25 CKl (HOCfc )=CKN 
DK=GKN 
CIFFI=DIFF 
IF ( IHI) GG TO 26 
IF (SLa^^) GQ TO 27 

CENTER WELL 

CKCEN2=CKCEN 

CKCEN=CKK 

JF1=JTPL 

JLCH=JTFL 

IF (JTPU.lt. JUPPER ) IHPCCE=0 
IF (KWELL.lt. NWELL) GC TO 28 
KWELL= 1 

IF ( JTPL.GT. JBOTT) L0NCCE=0 

ISFAP= 1 

RETURN 

UPPER WELL 



26 CFCLC=CKCEN 
JLPPER=JHI 
$FI=. FALSE. 
KWELL=2 
ISFAP=2 
RETURN 



r LOWER WELL 

27 CKCLD=CKCEN 
KWELL= 1 
$LCW=. FALSE. 

ISHAP=4 

RETURN 

CHECK FOR PODE IN LCWER DUCT 



28 KWELL=KWELL+1 
• ISFAP=3 

J1=JTPU 
CS1=DS 
J2=JTPL 
CALL WKS 

LCCIFF = CS-CFLCAT( LCMCCE + 1 I’pCPI 
IF (L0CIFF+1U.0DC>?^DTEST ) 29 ,30,3C 

29 KWELL=1 
JTPL=J2 
CS=CS1 
JTFU=J1 
RETURN 

SET UP FOR LOWER CHANNEL 

5C LCHCDE=LCMCDE+1 
jeCTT=JTFU 
FL=LODIFF 
4LCW=.TRUE. 

JTFU=J1 

CS=DS1 

JTFL=J2 
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RETURN 

C 

31 FCFR/lT (•0^'CDE ',14,' FCUNO BY CIPECT HI!.') 
ENC 
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SCEROUTINE SHAPED 

IMPLICIT ■REAL'i'O (A-H.V-Z), LCGICAL’?'^ ($) 




COMMON ELCCKS >>>>>> 



CCMiVQN 


/ 


D2 


/ 


ZZZ( IC2U , DVZ( LOCI) 


COMMON 


/ 


C4 


/ 


DRO,CRe,DRORb 


COMMON 


/ 


C5 


/ 


CMAX ,CMIN,CB, CESQ 


COMMON 


/ 


D6 


/ 


H,CF,hPLCT,DHhLF,0h2SIX,0H2C3, 


1 








0H2, hBOTT ,DH3 


COMMON 


/ 


07 


/ 


0K,CKM N ,DKMAX ,DKCLC,CKU,CKL, 


2 








DTEST,DSMAX, DE 


COMMON 


/ 


C9 


/ 


ZMAX,CSS,CECAY,CPVEL,CGVEL,GAM 


COMMON 


/ 


010/ 


WOCB ,NCCNP1 ,WOCM,I N,CVZB 


COMMON 


/ 


ow 


/ 


FQV ( 2G ) , FQ, W, N2 


COMMON 


/AWKB/ 


DKN,CINT,CIFF,CS, CM 


COMMON 


/ 


•I 1 


/ 


NUMV ,NMCD,NFREC,NMCCF,N,NPl , 


3 








NP21 ,NPTS,MODE, IT ,NP2 


COMMON 


/ 


12 


/ 


LOMCCE, IhMODE ,KWELL , IWELL ,NVvELL 


4 

COMMON 


/ 


13 


/ 


JHI , JLOW 

JUPPER,jeCTT, ISHAP, JTPL,JT PL 


CCMMON 


/ 


10 


/ 


NR£AD,NPRINT,NPUNCh 


COMMON 


/ 


L 1 


/ 


$GRAPh,$6CTPR,$CARCS,$CEL 


COMMON 


/ 


OTH/ 


DTHETU 




DATA >>>>>>»»>»» 
PROGRAM BEGINS >>>>> 



CSTART= I .CCO 

CE = (DK.VAX-CK) -MCKMAX + CK) 

IF (ISFAP,.EQ.3.ANC.IhiVCCE.GE.l) ISFAP=6 
IF (ISFAP..EQ.1.ANC.IFMQDE.GE.L) ISFAP=5 
$GC=. FALSE. 

I START = JTPL 
ISTlP=NP1 

GC TO ( l,l,l,2t3,3) , IShAP 

1 ITCP=1 
GC TO 4 

2 ITCP=JLCV« + l 
GC TO 4 

3 ITCF=JUFPER 

4 GC TO (7,5,6,7,7,61, IShAP 

5 ISTCP=JhI-l 
GC TO 7 

6 ISTCP=JECTT-1 

7 CONTINUE 

AVKZ = OI NT/ (CFLDAT( JTPL-JTPU ) 

ZJST=OS IN (OTHETU) *C START 

FACTOR=ZJST=^^( CVZ ( I START )-DVZ( ISTART-1) ) / ( A VKZ^4 . DO=^DH ) 

ZPJST = AVKZ*0CCS( DThETUI’^DST ART- FACTOR 

ZJ=ZJST 

2FJ=ZPJST 

FJ=DE-OVZ( ISTART) 

JSTART= ISTART 
JSTCF=I STOP 

CENOM=l .CDG-DH2SIX*(0E-DVZ( JSTART-in 
Z JM1=( ( I .0D0-0H2D3»FU )«ZJ-DH“ZPJ )/CENOM 



DEPTH LOOP BEGINS 

ZMAX=C.CD.O 
C 

DC 10 J=JSTART, JSTCP 
Z2Z( J)=ZJ 
FJ=DE-DVZ( J) 

ZJF 1= (2.000-DF2 + FJ)*ZJ-ZJMl 
ZuMl=ZJ 
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/!2J = CABS(ZJ) 

C 

IF (J-JTPU 9,9,8 
C 

e IF (AZJ.LT.ZFAX*CTEST ) GO TC 11 
IF (DAbS(ZJPl) .GT.AZ-J) GO TG II 
C 

G IF (AZ J .GT.OUPPER J GO TC 29 
IF (AZJ .GT .ZMAX) Zf^AX=AZJ 
1C ZJ=ZJP1 



DEPTH LOOP ENDS 



IF (ISTCP.EG.NPl) GO TO 16 

JSTART=JSTCP+1 

GC TO 12 

U JSTART= (J + jTPD/2 

12 2J=ZZZ( JSTART) 

DC 15 J=JSTART, JSTCP 
ZZ2( J)=ZJ 

CKZ=DSQRT{DA6S(DVZ(J)-DE) ) 
2J=ZJ^DEXP(-CKZ*CF ) 

IF (DAuS(ZJ) .LT.ZFAX-CTEST) GO TC 14 

13 CCNTINUE 

IF (ISTCP.EG.NPl ) GC TC 16 
J=JSTOP 

14 JSTART=J+1 
$GC = .TRLE.. 

DC 15 J=JSTAPT,NPTS 

15 ZZZ(J)=C.0D0 



CCMFUTE UPPER DECAYING SOLUTICN 



C 



c 

c 

c 



c 



16 IF (ISTART.EQ.ITOP ) GO TO 24 
ZJ=ZJST 
ZPJ=-ZPJST 
FJ=DE-0VZ ( ISTART) 

ZJM = ZJ-Dh’i'ZPJ-OH2=!'FJ*0 .SCO 
ISTAR1= ISTART+1 
JENC=ISTARl-ITOP 
ZRXDT = ZPAX=^'DTEST 



DC 18 J=1,JENC 
INC£X= 1ST ARl-J 
ZZZ( INDEX) =ZJ 
. FJ = DE-DVZ( INDEX) 

ZJF1 = (2 .0D0-Dh2’^FJ)*ZJ-ZJMl 
ZJM = ZJ 
AZJ=CAbS( ZJ) 

IF (AZJ-Zi'^XDT ) 20,20,17 

17 IF (0AdS(2;JP1)-AZJ ) 18,18,22 

18 ZJ=ZJP1 



19 IF (ITCP.EG.l) GO TO 24 
2C KST0F=INDEX-1 

DC 21 K=1,KSTCP 

21 ZZZ(K)=C.ODO 

GC TO 24 

22 JSTART= (J+n/2 
Zsl^ZZH ISTaRI-JSTART ) 



DC 23 J=JSTART, JENC 
INC£X=ISTAR1-J 
ZZZ( INDEX) =ZJ 

CKZ=DSGRT (CABS (DVZ( INDEX )-DE ) ) 
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Z J = ZJ*CEXP (-DKZ='<DH ) 

IF (CABS(ZJ) .LT .ZPXCT) GO 7C 20 
22 CCMINUE 
C 

GC TC 19 

NORMALIZE TFE FUNCTION 

2A CCCEF=l .ODO/ZMAX 

CC 25 K=l,NPl 
25 ZZZ(K )=ZZZ(K)*DCOEF 



COMPUTE TFE BOTTOM SCLUTICN 

IF ( .NCT.5B0TPR) RETURN 
EXFCN=DEXP (>CFB«DSCRT ( LVZb-CE ) ) 

22ZE0T=ZZZINP1)~DRCRE 

CC 26 J=NP2,NP2l 
ZZZ( J)=ZZZECT 

IF (DAbSIZZZBGT) .LT. 1. CD-60) GO TC 27 
26 ZZZEOT= ZZZECT=!=EXPGN 

RETURN 

C 

2? CC 28 I=J,NP21 

28 2Z2(I)=C.O 
C 

RETURN 

C 

c 

29 IF (DSTART .LT.l .CD-60) GO TO 30 
CSTART=CSTART=<=1 .OC- 10 

GC TO 7 

2C WRITE (NPRINT,31) 

J=1 

GC TO 27 



31 FORMAT Cl EXPONENT PRCTECTIOM EXCEDED IN SFAPE4',//, 
1 • RETURNED TC CALLING PROGRAM INCOMPLETE') 

ENC 
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SLERCUTINE OUTPUT (CUZ.NPTS) 

INPLICIT RE;:\L*8 (A-H, V-Z ) , LUGICAL^A ($) 

c 

CIPENSiCK DZZ(NPTS), DUZ(NPTS), CCC(NPTS) 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON ELCCKS >>>>>> 



COMMON 


/ 


D4 


/ 


DRO ,CRB,DRORB 


COMMON 


/ 


05 


/ 


CMiAX,CMIN,CB,CBSG 


COMMON 

1 


/ 


06 


/ 


H,CH ,HFLCT , DHHLF , DH2S IX , DH2 03 
DH2 ,HEOTT ,DHB 


COMMON 

2 


/ 


D7 


/ 


CK,OKMIN ,LKMAX ,DKCLC,CKU, DKL , 
DTEST,DSMAX,DE 


COMMON 


/ 


D9 


/ 


ZMAX , OSS, DECAY, DP VEL ,OGVEL ,GA 


COMMON 


/ 


OU 


/ 


FQV (20 , FC,W,W2 


COMMON 


/AWKB/ 


DKN ,CINT,CiFF ,DS,DM 


COMMON 


/HEAi 


D/ 


HED(20 


OCMMON 

3 


/ 


I 1 


/ 


NUMV ,NMGC ,NFREC, NMC0F,N,NP1 , 
NP21, NPTS, MODE, IT, NP2 


COMMON 


/ 


IQ 


/ 


NREAC, NPRINT, NPUNCH 


COMMON 


/ 


LI 


/ 


$GRAPH,$BCTPR ,$CARCS,$CEL 



C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< PROGRAM EEGINS >»» 
MODE OUTPUT SECTION 

1 WRITE (NPUPCH,6J PCC E , FC , CKK , CPV E L ,DGVE L , CS S , DECAY 
WRITE (NPUNCH,?) DLZ 

IF ( .NOT. SGRAPH) GC TO 2 
WRITE ( NPRINT, 9) PCCE 
RETURN 

2 WRITE (NPRINT, 10) 

RETURN 

SCUND PROFILE OUTPUT 

ENTRY PFOFIL (DZ2, COC, NPTS) 

WRITE (NPUNCh,3) HEC 

WRITE (NPUNCh,4) N PTS , CRO , CRB , C M I N , CB , H , HPL OT 
WRITE (NPLNCH,6> CZZ 
WRITE (NPUNCH, 8) CCC 
WRITE (NPRINT, 5) 

RETURN 
C 
C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< FORMAT STATEMENTS » 
C 

3 FORMAT (10A8) 

k FORMAT ( I1C,2F10.5,AF10.2) 

5 FORMAT (•CSCUND VELOCITY DETAILED PROFILE PLNCHED') 

6 format (I5,F5.0,G20.L2,2F10.2,2G15.8) 

7 FORMAT (1CF8.5) 

8 FORMAT (10F8.1) 

9 FORMAT CO OUTPUT FOR MODE', lA,' HAS BEEN PUNCHEC) 
IC FORMAT (•+', TL18, 'OUTPUT PUNCHED') 

END 
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SLEROUTINE INPUT2 

1^FLICIT REAL”8 LGGICAL^^ ($) 



LAST CHARGED 3 AUGbST L973 



COMMON BLOCKS >>»» 



CCMPGN 


/ 


C3 


/ 


CC(5C),CEP(5p) 


CCMMCN 


/ 


04 


/ 


DRC ,LRB,CRORB 


CCMMON 


/ 


C5 


/ 


C.MAX ,CMIN,C6,CBSQ 


CCMPON 

1 


/ 


06 


/ 


H,CH,HFLCT,DHHLF,CH2SIX,DH2C3 
DH2 ,HBOTT ,DHB 


CCMMON 

2 


/ 


D7 


/ 


DK,DKMIN ,CKMAX,DKGLC ,DKU,DKL , 
DTEST,OSMAX,D£ 


CCMMON 


/ 


DW 


/ 


FQV(2C) ,FG,U,1'2 


CCNMCN 


/H 


!EAD/ 


HEC(20) 


CCMMON 

3 


/ 


I 1 


/ 


NUMV,KMOD,NFREG,NMGCF,N,NP 1 , 
NP21 ,NPTS ,MODE, IT ,(\F2 


CCMMON 


/ 


IC 


/ 


NREAD,NPRINT,NPUNCH 


CCMMON 


/ 


LI 


/ 


$GRAPH, $BGTPR,$CABCS ,$CEL 



CATA >>>>»>»>>»» 



CATA CCNVPT / 1 . 000 C COOCCOCOOOOO CO/ 



PROGRAM BEGINS >>>» 



^ REAO IN THE HEADING 

AND RUN PARAMETERS 

READ (NREAC,9J HEO 
RRITE (NPRINT,13; HED 

READ (NREADtlO) NUPVt NMCD 1 1 EX ,N , ICLT , !DEV 

— . — : test the run parameters 

IF (N.GT.IOOO) N=1000 
IF (NUMV.GT.5C) NLMV=50 
IF (NMQC .GT . 1001 NMCD= 100 
IEX=IAeS(IEX) 

IF (IEX.GT.14) IEX=14 
DTEST = LC.ODO-^v(-IEX) 

IF (I0UT.LE.2) $GBAPH=.TRUE . 

IF (( I0UT.E0.2) .OR. nCUT.EQ.5)) $CEL=.TRUE. 

IF (CIOUT.EQ.l) .OR.( I0UT.EQ.4) ) $CAR0S= . TRLE . 

IF (lOEV.NE.O) NPUNCH.= iCEV 

NP1=N+1 

NP21=N+21 

READ (NREADtll) H , CB , DPC , DRB ,HP L CT 

read IN THE FREQUENCIES 

READ (NRcAD.lO) NFPEG 

IF (NFREC .GT.20 ) NFPEC = 20 

READ (NRcAOtll) ( FC V ( I ) ,I =l ,NFR EQ ) 

FC = FQV( 1) 

CONVERT AND DISPLAY 

THE RUN PARAMETERS 

IF (NUMV) 1,2,2 

1 CCNVRT=C .3C48 
NLMV=-KLPV 

H = H=?‘CCNVRT 

CE=Co*CCNVRT 

HPLCT=HPLCT^CONVRT 

2 V»PITE (NFRINT,L4) 

WRITE (NPRINT,15) ( FQV (I ) , I = 1 , NFR EG ) 

WRITE (NPRINT,16) NMQD,IEX 
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WRITE (WPRINT,17) H.HPLCT 
WRITE (RPR1KT,18) C0 ,N 
WRITE (WPRINT,19J 
IF UGRAPh) WRITE (\PRINT,22) 

IF (5CARCS) WRITE (i\FPlKT,21J 
IF UCEL; WRITE (NPRINT,20) IDEV 
WRITE (NPRINT,23) FED 

, jn TEE DEPTHS 

A^D SOUND VELOCITIES 



CC 2 

READ (NREACtlU CEPI.CCI 
DEP I = DEPI*CGNVRT 
CC I=CCI’i'CONVRT 
CC ( I ) =CCI 
CEPI I)=DGPI 

WRITE (NPRINT,12) DEPI,CCI 
3 CCNTINUE 



C 

C 



CF = F/CFLCAT(N ) 

CRCR6=LhC/CRE 
CF2 = DH--CH 
CESC = CB=?CB 

DF2SIX = Cl-2-^^0. 1666c66t666667 
CF2D2 = DH2S I)(i'2.CDC 
CFFLF=t.E=’^C .500 
IF lOEPI-H) 4,6,5 
A NLR'V=NUPV+1 

5 CC (NUMV )=(E-DEPI 1=^0 .016 500 + CCI 
CEP(N'UMV)=H 

6 NFTS=NP2L 

IF (HPLCT-h) 7,24,8 

7 FPLGT=H 

24 iECTPR= .FALSE. 

NFTS=NPl 

RETURN 

e HECTT=hPLG7-H 
CFE=HB0TT^-0.50-l 
RETURN 



•CGRPUTE SCRE CONSTANTS 



FORMAT ST/TEMENTS >> 



s 


FORMAT 


(10A3) 








IC 


FORMAT 


(5110, 


12) 






1 1 


FORMAT 


(8F10. 


3) 






12 


FORMAT 


( 'O' , 


T25, 


F6.1, T57, F6.1J 




13 


FORMAT 


Cl', 


10A8) 






14 


FORMAT 


co». 


T34, 


'RUN PARAMETERS ; ' ) 




15 


FORMAT 


CO', 


T2 3, 


'FREGUENCIES :',(T42, 


F6 . 1 , ' HZ' 


U 


FORMAT 


CO', 


T24, 


'MODES REQUESTED : ' , 


14 , // , 



/n 



1? 

le 

IS 

20 

21 

22 

23 



1 T24 

FORMAT 
I T22 

FORMAT 

« 



,• ITERATION LIMIT : 10*« 
(•C* , T27, 'BOTTOM DEPTH 
, .'MAX DEPTH PLOTTED 



( -' 
• « , 



f 12, 

F3. 1, 



) ' ) 

METERS 



Ffi.l, ' METERS' 



CO' 



T18 



1 ' METERS/SEC ,// , 
FORMAT ('O', T22, 

1 'NUMBER' ,//,T42 , • 



FORMAT CO 
FORMAT CO 
FORMAT ( ' ' 

FORMAT C C 
I 'DEPTH, 
ENC 



T42, 

T42, 

742, 

luAB 



'BOTTOM SOUND VELOCITY :',F8 
Ti7,'FUMSER OF GRID POINTS : 
'OUTPUTS RcwUESTED : H0RI2 WAVE 
GROUP AND PHASE VELOCITIES',//) 
'OUTPUT ON FILE FT ' , I 2 , ' F CC 1 ' ) 

• FUNCHEC CARC OUTPUT ' } 

•SOLNC VELOCITY AND MODE 
//, 733, 'VEL 



) 

I , 

15 ) 

I 



/ / , 



GRAPHS ' ) 

;CITY PROFILE ' ,// ,T22, 
METERS', T44, 'SOUNC VELOCITY, M/SEC, /) 
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SLBROUTIKe MOOPLT 

IMPLICIT REAL-^8 (A-H, V-Z ) , L0GICAL1=4 {$) 
CIPENSICN RANGE (4) 



4 

=4^ THIS SUBROUTINE PLOTS THE MCCE SHAPE CN THE PRINTE 
^ UTILIZING THE STANCARO (AT NFS) ROUTINE UTPLCT 



CCNPON / Cl / 
CCPMON / 02 / 
CCPMON / C5 / 
CCNPON / 06 / 

i 

‘CCPMON / D1 / 
) 

'CCPMON / 03 / 
CCPPQN / 09 / 
CCKPON/OKMAP/ 
CCPMGN / OW / 
CCPMON /AV'JKB/ 
CCPP.ON /HEAD/ 
CCPMON / II / 
5 

COMMON / 10 / 



VELCdOZl), DEPTH (1021) 

ZZZ( 1C21) , DVZ( LGC IJ 

CMAX ,CM IN ,Cd, CBSQ 

H ,Ch jt-FLCT ,DHhLF ,CH2SIX, DH2C3 , 

0H2iHB0TT,DHB 

OK ,CKMIN ,CKMAX ,DKCLC ,CKU, OKL , 
DTEST ,DSMAX,DE 
DIFF I 

ZMAXjCSS, CECAY ,DPVEL,OGVEL ,GAN 
OKI ( ICO) 

FQV(20) ,FQ,W,U2 
DKN',DINT,CIFF,CS,CH 
HED( 20 

NUMV,NM0C,NFREC,NMCCF,N,NP1, 
NP21,NPTS jMOOE, IT ,NP2 
NREAutNPRiNT » NPUNlH 



WRITE (NPRINT,!) HEC,FC 
WRITE (NPRINT, 2) MODE, OK 
WRITE (NPRINT, 3) 

CALL UTPLGT ( ZZZ , C E P TH , NPT S , RANG E , 2 , 0 , P BC T T ) 
WRITE (NPRINT, 4) CPVEL , CSS , DEC AY , CGVEL , I T , C IFF I 
RETURN 



ENTRY CZFLCT 
f BCTT=5NGL(H) 

WRITE (NPRINT, 5) HEC 
WRITE (NPRINT, 6) 

WRITE (NPRINT, 7) 

RANGE! 1 ) = SNGL(Cd) 

RANGE(2)=SNGL(CMIN) 

RANGE(4 ) = 0 .0 
RANGE(3)=SNGL(HPLCT) 

CALL UTPLCT ( VELC , CEPTH ,NPTS, RANGE ,2 ,0 , PECTT ) 

RANGE! I )=1 .0 

RANGE(2)=-l-.0 

RETURN 



1 FORMAT (M', 10A8, ' FREQUENCY : ', F6 . 1 , • HZ') 

2 FORMAT ('O', T25, 'MODE', 14, • : K =*, F16.12) 

3 FORMAT ('O', // , T23, 'NORMALIZEC £ IGEN FUNCT ICN ' , 

I • VS CEPTH •) 

4 FORMAT COPHASE VELCCITY F7.1, • M/SEC. OSS', 

4^ I !*,G15«7, 

1 4X, ' BOTTOM LEAKAGE COEF : • , G15.7, //, 

$ • GROUP VELOCITY 

2 F7.1, • M/SEC. ITERATICNS 15, ICX, 

3 'DIFFERENCE FUNCTION G16.E ) 

£ FORMAT Cl', 10A8 ) 

6 FORMAT CC, //, T24, 'GRAPH OF SCLND VELOCITY', 

I ' VS DEPTH' ) 

7 FCPMAT CC, T19, 'CEPTH IN METERS.', 

1 ' SCUND VELCCITY IN METERS/SEC.') 

END 
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SLBROUTIKE LIKEAR 

IMPLICIT R£AL*8 (A-h, V-Z ) t L0G1CAL^=A ($) 



C 



CCMMON ELCCKS >>>>>> 



CCRRCM 

CCMRON 

CCPKGN 

CCRNON 

CCRRON 

1 

CCMMGN 

CCNMON 

CCPHON 

■a 

■'CCRMON 



/ Cl / VELCdUZl), CEPTH (1021) 

/ 02 / ZZZ( 1C21) , DVZ(lOOl) 

/ 03 / CC(5C),DtP(50 ) 

/ 05 / C«AX,CPIN,CS,CeSG 
/ 06 / H,Ch,hPLCT,0HHLF,0h2SIX,0H2C3, 
Dh2,hEGTT,ChS 

/ 05 / ZMAX,CSS,CECAV ,OPV£L,CGVEL,GAR 
/ DW / FQV(20) ,FQ,W,Vh 2 
/ II / NUMV,NMCC,NFRE0,IMRCCF,N,NP1, 
NP21 ,KPTS,MC0E, IT ,RF2 
/ LI / $GRAPF, $BCTPk , $CAR0S ,1CEL 



C 



PRQGRAP BEGINS >»» 



CPAX=C. COO 
CP IN=CB 



CL = C .OCO 
CL = DEP( 2 ) 

VL = CC( I ) 

VL = CC (2 ) 

CLPCU=DL-OU 

VLPVU=VL-VU 



1=2 



CEPTH LGCP BEGINS 



CC ^ J=1,NP1 
CEPJ=OH*CFLCAT( J-1 ) 

0£PTH( J ) =DEPJ 
IF (DEPJ-DL } 2,2,1 
1 1 = 1+1 
Cl = 0L 
VL=VL 
CL = CEP( I ) 

VL = CC( I ) 

CLP0U=0L-0U 

VLPVU=VL-VU 

2 VELOJ = ( OEPJ-OU)’^VLPVL/OLMOU+VU 

IF (VELCJ.GT.CPAX) CMAX=VELCJ 
IF (VELCJ-CMIN) 3,4,4 

3 CMN=VELCJ 

4 VEL0( J)=VELOJ 



CEPTH LGCP ENOS 

IF ( .NCT.$BGTFR) GG TO 6 
NF2=NP1+1 

ASSIGN 6CTTCM VELCCITY VALUES 



CC 5 J^1,2C 
VELGINP 1 + J)=CB 

CEPTH(NP1.+ JJ = F+0HB*CFLCAT( J-1) 
5 CCNTINUE 



6 RETURN 
ENC 
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c 

c 

c 

c 

c 



c 

c' 



c 

c 

c 

c 



c 

c 

c 

c 



c 

c 

c 



SLfcRCUTINE WKB 

IMPLICIT REAL^e ( A-htV-Z ) ,LOGICAL^'i ( S) 



LAST CHANGED 2 7 SEPT 73 



AIRY PHASE B.C 



COMMON BLCCKS >»>» 



CCPNCN / C2 / 
CCMMON / DA / 
CCMNON / C6 / 

^CCMHON / 07 / 
) 

'CCMMON / 010/ 
CCMMQN /AWKB/ 
CCMMQN / II / 

X 

'CCRMQN / 12 / 



COMMON 

COMMON 



13 / 
LI / 



COMMON / PI / 
COMMON /^DTH/ 
DATA OEXRU/l. 



ZZ2(1C21I, CVZdOOlJ 
ORO ,CRB,DRORB 

H, CH ,HPLCT,DHHLF, DI-2SIXtDH2C2, 
0H2 ,HEGTT ,DHB 

OK fOKMIN tOKMAX ,0 KC LC , OKU » OK L , 

DTEST ,CSMAX,CE 

WOCB ,V\CCNP1 tWCCMI N tCVZB 

DKK,DINT,DIFF,CS,CM 

NUMV ,NMCC ,NFREC,NMCCF,N,NP 1 , 

NP21 ,NPTS .MODE, n .KF2 

LOMCDE, IHMODE.KWELL.IWELL .^l^ELL, 

JHl , JLQK 

JUMPER, J60TT, liHAP .JTPL.JT PL) 
SGRAPF, $60TPR,$CARLS,S.CEL 
OP I, CPI 02 ,0PIC4,DFIBT2,D2PI 
DTHETL 

702/ .CWELL/6 .0/ 



PROGRAM EEGINS »>» 



CINT=0.C0C 
5TCP=. FALSE. 

JTFL=1 

C E= ( DKM AX-DKN )--!={ DKMAX + CKN) 
CKE=DSC;RT (CVZB-CE ) 

IWELL=1 
ISIG=-l 
CARG=C. CCO 
DAPGP=C.CDC 



CC 21 J=1,NP1 
CKZ2=DE-DVZ( J ) 



•CEPTH LCCP BEGINS 



•CHECK FOR CONDITIONS 
AND TYPE CF SOLUTION 



IF (DKZ2) 1,2,3 
1 IF (ISIG) 4,5,6 
IF (ISIG) 7,21,S 
IF (ISIG) 10,19,20 



IMAGINARY REGION 



CARGP=DSCRT(-DKZ2) 
CARG=CARG.+ CARGP 
GC TO 21 



5 ISIG=-1 
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GC TO ^ 

C 

e C/if<GP = CSCRT(-CK22 ) 

CX=DKZ/ (CKZ+OARGP ) 

CIM=LIM + DKZ^(DX-1 .0C0)’f'0.5D0*Ch 
ISIG=-l 
JTFL=J-i 

CZRG=(2.0C0-0XJ^0ARGP«C.5D0 
GC TO 21 

2ER0 VERTICAL VvAVE NUMBER 



e ISIG=0 
GC TO 21 

S IEACK=1 
GC TO 8 

real REGICN 

1C ISIG=1 

IF (ITQPJ GO TO 11 
JTFU=J 

11 CKZ=DSGRT{CKZ2) 

CX=CARGP/ (CKZ+OARGP ) 

C/!RG=CARG^CARGPV(DX-1 .OCO )*0.5DC 
IF ($TOPi GO TO 13 

12 C2=2.CCC-D£XP(2. ODC-CARG+DH ) 

CTFETU=CATAN( (D2-1 . COO ) / ( 02+ 1 .0 CC )) 
CINT=OTHETL 

$TCP=.TRUE. 

GC TO 18 

13 CT 1=DMCC(0INT,02PI ) 

IF (OARGvCh-CWELL ) 15,1A,1A 

lA IF (IWELL.EC.KWELL) GC TO 28 
CARG=C. CCC 
IRELL=IRELL+1 
JTFU=J 

DThETU=CPIOA 
CINT=OTFETU 
GC TO le 
15 A=CSIK(CT1) 

E=CCGS(CT1) 

CEL=C£XF(OARG*OH) 

CEPL=C£XP(-DARG-i'CI- ) 

AC = DEL-= ( A+E ) + C£ML- ( A-B ) 

EC = 0EL'MA + B)-0EML'"'( A-E) 
CINTP=DINT-0T1+CATAN2( AG.BO J 
IE IF (OINTP.GE.OiNT-CTEST) GC TO 17 
• CINTP=D INTP+0PI02 
GC TO 16 
1? CINT=OINTP 

1£ C IKT=0INT + 12.000-CX)*0KZ*C.500*DH 
CARG=0.CC0 
GC TO 21 

19 CKZ=OSURT(DKZ2) 

1SIG=+1 

C INT=OINT.+ OKZ=»OH 
IF ($TCF) GO TO 13 
JTPU=J 
GC TO 12 

2C CKZ=CSQRT(OKZ2) 

C IRT = 0I NT + CKZ*0H 

21 CONTINUE 



DEPTH LCCF ENDS 



15A 



ooo o ooo non 



IF (DKZ2) 22,27,25 

DECAYING SCLCTICN 

22 CKS=DARGF 

CAFG= ( CAP,G-C.5D0*DARGP )=?DH«2.0D0 
IF (CARG-DEXPU) 2b,24,2A 

22 D2=(DKS-DRCRb»'DK8 ) ’f^OEXP ( -DARG )-^2 .0 CO/ ( DKS + C PCR E=pOK B I 
DThETL=CATAN( ( l.0D0+D3)/( 1.000-D2) ) 

GC TO 27 
Ih CTFETL=DPI04 
GC TO 27 

ECTTOM REFLECTED 

25 JTFL=NF1 

IF (DKB.EQ.O.ODO) GC TO 26 
CTFETL=DATAN(DKZ/ (CK8*DRORBn 
GC TO 27 

26 CTFETL=CFIC2 

— COMPUTE DIFFERENCE 

27 CS=CINT+CTFETL 

C INT=DINT-DThETU 
CIFF = DS-DM->DPI 
NVnELL=I V\ELL 
RETURN 

2€ UELL = I VsEU+L 
GC TO 2A 

END 
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